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CHAPTER  1 


INTRODUCTION 

Background 

The  growth  of  many  cities  in  the  United  States  has  resulted  in  urban 
development  in  many  natural  floodplains.  This  development  often  increases  runoff 
from  rainfall  which  in  turn  increases  river  flow  rates.  These  increased  volumetric  flow 
rates  have  resulted  in  floods  which  have  caused  millions  of  dollars  in  property 
damages  and  loss  of  life. 

Engineers  are  often  asked  to  propose  methods  to  reduce  flood  damages  and 
prevent  loss  of  life  on  developed  floodplains.  Construction  of  flood  control  channels 
is  a  popular  method  to  increase  the  capacity  of  the  natural  channel  and  reduce 
flooding.  These  channels  are  typically  concrete  lined  or  constructed  as  composite 
channels  with  the  invert  concrete  and  the  side  slopes  grouted  stone.  If  the  natural 
topographic  relief  is  sufficient,  these  lined  channels  usually  are  designed  with 
hydraulically  steep  slopes  to  convey  supercritical  flow  over  most  of  the  improved 
reach.  Froude  numbers  in  these  man-made  channels  commonly  range  from  1.2  to  2.0. 

Hydraulic  engineers  many  times  use  the  term  "high-velocity  channels"  when 
referring  to  a  lined  flood  control  channel  designed  to  discharge  supercritical  flow.  The 
designer  of  these  high-velocity  channels  in  developed  areas  is  faced  with  many 
questions  which  are  not  easily  answered.  The  main  concern  is  the  depth  of  flow  in  the 
channel  for  the  design  discharge.  The  depth  must  be  known  to  determine  sidewall 
heights  and  minimum  bridge  soffit  elevations.  Determining  the  depth  of  flow  is 
complicated  by  side  inflows  and  boundary  features  such  as  contractions,  expansions. 
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curves,  and  obstructions  to  the  flow  such  as  bridge  piers,  and  vehicle  access  ramps. 
These  boundary  features  in  a  supercritical  channel  cause  flow  disturbances  which  can 
result  in  a  significant  increase  in  the  local  flow  depth.  Figure  1.1  shows  the  flow 
conditions  in  a  high-velocity  channel.  The  oblique  standing  waves  in  the  foreground 
are  generated  at  the  tails  of  piers  located  just  upstream  of  the  reach  shown  in  the 
photograph.  An  accurate  prediction  of  the  water  surface  shape  (i.e.  variations  in  local 
depth)  is  essential  in  the  successful  design  of  a  high-velocity  channel. 
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The  cross  sectional  shape  of  high-velocity  channels  is  generally  either 
rectangular  or  trapezoidal,  the  choice  of  which  is  dictated  by  optimization  of  hydraulic 
and  economic  concerns.  The  economic  tradeoff  of  rectangular  versus  trapezoidal 
channels  is  that  rectangular  channels  require  less  real  estate;  however,  construction  of 
rectangular  channels  is  more  costly  than  trapezoidal  channels.  The  structural  cost  of 
vertical  sidewalls  is  often  greater  than  the  cost  of  acquiring  the  additional  real  estate 
required  to  convey  the  flood  flow  in  a  trapezoidal  channel.  Trapezoidal  channels  can 
be  constructed  by  using  levees  on  both  sides  of  the  channel  and  then  paving  the 
channel  side  slopes  and  invert.  Channel  cross  section  shapes  are  selected  based  on 
adequate  hydraulic  capacity  in  conjunction  with  the  minimum  cost  of  construction  and 
maintenance.  The  total  cost  includes  costs  of  right  of  way,  the  cost  of  construction 
and  modification  of  structures  such  as  bridges,  and  maintenance  costs  associated  with 
the  removal  of  sediment  and/or  debris. 

The  flow  in  trapezoidal  high-velocity  channels  differs  markedly  from  flow  in 
rectangular  channels.  Straight  reaches  of  rectangular  channels  produce  a  fairly 
uniform  distribution  of  velocity  at  a  cross  section,  although  the  actual  distribution  is 
influenced  by  the  channel  boundary  and  free  surface  drag  and  therefore  depends  on  the 
channel’s  width-to-depth  aspect  ratio  (Rouse  1961).  However,  the  velocity  variations 
across  a  section  of  straight  trapezoidal  channel  are  large.  Flow  along  the  side  slopes 
near  the  water’s  edge  is  severely  retarded  by  bed  drag  as  the  flow  depth  decreases. 

The  local  Froude  number  near  the  side  boundary  approaches  zero.  The  local  flow 
regime  at  a  cross  section  of  a  trapezoidal  high-velocity  channel  can  be  supercritical  at 
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the  center  and  subcritical  near  the  water’s  edge.  The  flow  over  a  portion  of  the  side 
slope  may  be  at  or  near  critical.  This  mixture  of  subcritical  and  supercritical  flow 
regimes  can  lead  to  flow  instabilities  and  associated  fluctuations  in  depth  and  velocity. 

Hydraulic  Design  of  High-Velocity  Channels 
Supercritical  flow  in  channels  is  characterized  by  standing  waves  created  by 
even  small  changes  in  the  sidewall  alignment  (figure  1.2a).  Oblique  standing  waves 
are  created  at  the  beginning  of  channel  bends  (PC,  figure  1.2b)  where  the  depth  along 
the  outer  wall  increases  as  the  wall  exerts  a  force  on  the  fluid.  This  depth  increase  is 
also  realized  at  the  beginning  of  a  channel  contraction.  The  flow  along  the  inside  wall 
at  the  point  of  curvature,  PC,  responds  as  at  a  channel  expansion  resulting  in  a  depth 
reduction.  However,  away  from  the  walls  the  flow  remains  unchanged  until  the  forces 
generated  at  the  boundaries  are  present  in  the  form  of  pressure  differences.  The 
standing  wave  pattern  downstream  of  the  initial  wave  intersection  point  is  not  a 
straight  line,  because  the  flow  at  this  point  is  moving  along  a  curved  path.  Reflections 
of  these  waves  continue  over  long  distances  downstream  of  the  curve.  The  water 
surface  setup  or  superelevation  in  rectangular  channel  bends  (figure  1.3  a)  is 
significantly  greater  in  supercritical  flow  than  in  subcritical  flow  and  the 
superelevation  in  trapezoidal  channels  results  in  runup  on  the  side  slope  at  the  outside 
of  the  bend  and  drawdown  on  the  inside  side  slope  (figure  1.3b). 

Before  the  classic  papers  composing  the  1951  symposium  on  high-velocity 
flow  were  published  by  the  American  Society  of  Civil  Engineers  (Ippen  1951,  Knapp 
1951,  Ippen  and  Dawson  1951,  and  Rouse  et  al.  1951)  each  component  of  a 
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supercritical  channel  design  was  tested  in  a  physical  hydraulic  model.  The  symposium 
papers  provided  fundamental  descriptions  of  the  mechanics  of  supercritical  flow  in 
man-made  channels.  These  descriptions  provided  design  guidance  for  engineers  who 
designed  high-velocity  channels. 


FLOW^ 


STANDING  WAVES 


a.  Rectangular  channel  contraction 


b.  Rectangular  channel  bend 

Figure  1.2.  Schematic  plan  views  of  standing  wave  patterns  in  supercritical  flow. 
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a.  Rectangular  channel 


b.  Trapezoidal  channel 

Figure  1.3.  Cross-sections  showing  superelevation  in  a  channel  bend  (r  is  the  inside 
bend  radius). 

Ippen  (1951)  details  the  treatment  of  oblique  standing  waves  for  supercritical 
flow  at  small  changes  in  the  vertical  sidewall  boundary  alignment.  Ippen  gives  the 
oblique  wave  angle  and  the  resulting  wave  height.  Water  surface  disturbances  are 
analyzed  by  either  assuming  constant  specific  energy  or  by  considering  the  energy  loss 
across  the  disturbance.  The  first  assumption  is  appropriate  for  gradual  water  surface 
changes  while  analyses  of  large  standing  wave  fronts  must  consider  energy  losses.  In 
his  analysis  Ippen  makes  the  assumptions  that  the  flow  is  steady,  the  pressure 
distribution  is  hydrostatic,  the  velocity  does  not  vary  in  the  vertical,  and  the  channel  is 
frictionless  or  that  the  friction  slope  is  parallel  to  the  channel  slope. 
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Ippen  and  Dawson  (1951)  demonstrate  the  applicability  of  the  findings  of  Ippen 
(1951)  to  design  problems  of  channel  contractions.  This  is  accomplished  by 
comparing  analytical  solutions  to  laboratory  data.  These  comparisons  reveal  that  phase 
errors  exist  between  the  theoretical  and  observed  depth  profiles  along  a  contraction. 
The  theoretical  waves  are  located  upstream  of  the  observed  waves.  Ippen  and  Dawson 
attribute  this  discrepancy  to  the  hydrostatic  pressure  distribution  assumption  used  in 
the  analytical  model.  However,  they  note  that  the  analytical  solution  is  generally 
adequate  for  chaimel  design  purposes.  They  also  provide  design  guidance  for  ways  to 
reduce  waves  downstream  of  a  contraction. 

Knapp  (1951)  gives  a  general  description  of  supercritical  flow  in  rectangular 
and  trapezoidal  channels.  Knapp  provides  a  means  of  estimating  the  superelevation 
and  standing  waves  within  a  curve  and  points  out  that  waves  generated  in  a  curve  are 
continued  downstream.  Knapp’s  assumptions  include  those  given  by  Ippen  with  the 
addition  that  the  velocity  is  constant  across  the  cross-section.  Knapp  states  that  in 
trapezoidal  channels  the  depth  is  not  constant  and  that  the  wave  celerity  therefore 
varies  throughout  the  cross  section  making  his  findings  much  less  applicable  to 
trapezoidal  channels.  He  cautions  channel  designers  that  the  additive  effects  of  curves 
separated  by  short  distances  are  not  treated  in  his  analyses  since  his  method  neglects 
the  nonuniform  velocity  distribution  across  the  beginning  section  of  the  downstream 
curve.  The  skewed  velocity  distribution  entering  the  downstream  curve  is  generated 
by  flow  through  the  upstream  curve. 

Considerations  are  given  by  Rouse  et  al.  (1951)  for  the  design  of  channel 
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expansions  in  the  presence  of  supercritical  flow.  Their  assumptions  include 
hydrostatic  pressure  distribution,  negligible  friction  losses,  a  flat  channel  (no  slope), 
and  vertical  sidewalls.  Rouse  et  al.  provide  graphical  means  for  estimating  the  surface 
configuration  of  the  negative  wave  generated  at  abrupt  expansions  and  provide 
guidance  for  the  design  of  efficient  channel  expansions  which  eliminate  disturbances  at 
the  downstream  end  of  the  transition. 

One-Dimensional  Analysis 

Typically,  in  engineering  applications,  open  channel  flow  is  analyzed  using 
empirical  one-dimensional  hydraulic  equations  where  the  velocity  and  the  water 
surface  elevation  are  assumed  constant  over  a  cross  section.  Fluid  frictional  losses  in 
uniform  flow  in  straight  prismatic  channel  sections  are  represented  with  either  the 
Manning  or  Chezy  equation.  Empirical  expressions  have  been  developed  to  estimate 
the  head  losses  (Henderson  1966)  and  superelevation  of  the  water  surface  in  channel 
curves  (Woodward  and  Posey  1941).  Flow  obstructions  such  as  bridge  piers  have 
been  studied  and  parameters  such  as  pier  loss  coefficients  (USAE  Los  Angeles  District 
1939)  have  been  quantified  for  simple  geometric  configurations.  Hydraulic  research 
has  resulted  in  design  criteria  for  channel  width  transitions  (Brater  and  King  1976,  US 
Bureau  of  Reclamation  1967,  Ippen  and  Dawson  1951).  The  general  location  of  a 
hydraulic  jump  in  a  straight  prismatic  rectangular  channel  is  determined  by  backwater 
computations  from  a  downstream  control  point  and  forward  step  computations  from  an 
upstream  control  point  until  the  momentum  equation  of  the  hydraulic  jump  is  satisfied. 

A  more  complete  analysis  must  be  conducted  for  complex  geometries  having 
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variable  cross  sections  or  obstacles  in  the  flow.  This  often  involves  using  a  physical 
hydraulic  model.  The  time  and  cost  requirements  of  physical  models  are  necessitating 
the  continued  development  and  use  of  numerical  flow  models  for  analyzing  flow 
patterns  in  complex  channels. 

Two-Dimensional  Analysis 

One-dimensional  analyses  are  not  adequate  to  describe  flow  patterns  in  which 
the  depth-averaged  velocity  varies  across  the  cross  section,  particularly  in  the  analysis 
of  supercritical  flow.  Moreover,  the  one-dimensional  analysis  fails  to  describe  the 
additive  effect  of  disturbances  created  at  various  points  along  a  channel  reach  such  as 
two  curves  separated  by  a  short  distance.  One-dimensional  modeling  of  channel 
features  such  as  bridge  piers  and  transitions  requires  knowledge  of  applicable  loss 
coefficients.  While  significant  research  has  been  directed  toward  quantifying  these 
coefficients,  general  energy  loss  relations  have  not  been  documented  for  complicated 
geometries  such  as  at  bridges  having  complex  pier  configurations,  at  skewed  bridges, 
or  at  bridges  crossing  channel  bends.  A  more  complete  description  of  the  flow  is 
provided  by  a  depth-averaged  two-dimensional  analysis  where  the  horizontal  velocities 
and  the  water  surface  elevation  vary  in  both  the  longitudinal  and  transverse  directions. 
Two-dimensional  modeling  of  flow  around  bridge  piers  and  through  transitions  does 
not  require  empirical  loss  coefficients  but  rather  only  requires  adequate  boundary 
descriptions. 

Depth-averaged  two-dimensional  flow  analysis  is  conducted  using  the  shallow 
water  equations.  The  shallow  water  equations,  sometimes  referred  to  as  the  de  Saint- 
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Venant  equations,  have  been  used  to  describe  two-dimensional  steady  and  transient 
free-suiface  flows.  The  principle  assumption  in  this  system  of  non-linear  hyperbolic 
partial  differential  equations  is  that  vertical  accelerations  are  negligible,  which  is 
equivalent  to  assuming  that  the  pressure  distribution  is  hydrostatic.  This  assumption  is 
adequate  when  the  water  surface  curvature  is  small  compared  to  the  depth  of  channel 
flow.  The  channel  slope  is  often  considered  geometrically  mild  such  that  sin  a  ~ 
tan  a  »  a,  where  a  is  the  channel  slope.  This  assumption  is  excellent  for  a  <  5 
degrees  and  so  is  applicable  to  man-made  channels  passing  through  alluvial  valleys 
having  maximum  slopes  of  about  0.015  («  0.9  degrees).  This  geometrically  mild  slope 
assumption  which  complements  the  hydrostatic  assumption  is  distinguished  from  a 
hydraulically  mild  slope  which  conveys  uniform  flow  in  the  subcritical  regime.  That 
is,  the  shallow  water  equations  assume  the  slope  is  geometrically  mild  although  it  may 
be  hydraulically  steep.  The  shallow  water  equations  are  capable  of  modeling  both 
subcritical  and  supercritical  flows.  The  depth-averaged  two-dimensional  flow 
equations  used  in  this  dissertation  are  presented  in  Chapter  3. 

Analytical  solutions  of  these  equations  are  not  available  for  two-dimensional 
problems  unless  simplifying  assumptions  are  made.  Generally,  steady  fluid  flow  is 
assumed  and  often  the  depth  and  cross  sectional  velocity  distributions  are  prescribed  as 
boundary  conditions.  Because  there  are  no  general  analytical  solutions  available  for 
the  shallow  water  equations  for  general  boundary  geometry  and  transient  channel 
inputs,  the  equations  are  solved  numerically.  Until  recently,  the  majority  of  fluid 
dynamics  and  hydraulic  problems  have  been  numerically  modeled  using  the  finite 
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difference  method.  In  particular  the  MacCormack,  Lax,  Gabutti,  and  Beam  and 
Warming  schemes  have  been  applied  to  the  shallow  water  equations  in  modeling 
supercritical  flow  in  rectangular  channels  in  which  fluid  flow  shocks  are  created 
(Jimenez  and  Chaudhry  1988;  Fennema  and  Chaudhry  1989,  1990;  Elliot  and 
Chaudhry  1992;  Molls  1992;  and  Younus  and  Chaudhry  1994).  These  schemes  are 
explained  in  texts  such  as  Anderson  et  al.  (1984).  The  finite  element  method  is 
employed  because  of  its  ability  to  describe  the  boundaries  of  complex  geometries 
accurately  and  its  straight  forward  manner  of  satisfying  boundary  conditions 
(Katopodes  1984a  and  1984b,  Olsen  1974). 

In  the  solution  of  two-dimensional  open  channel  flow  equations  the  local  cross- 
section  and  depths  are  unknown  quantities.  Numerical  solution  of  the  flow  variables 
treats  the  flow  field  as  a  number  of  discrete  flow  regions.  If  the  channel  sidewalls  are 
vertical,  the  plan  view  domain  is  well  defined  such  that  one  can  easily  discretize  the 
domain.  A  computational  grid  can  be  constructed  to  represent  the  channel.  However, 
if  the  channel  has  sloping  sidewalls,  the  plan  view  of  the  flow  domain,  as  delineated 
by  the  water  surface/bank  interface  (i.e.  the  "waterline"),  is  not  known  a  priori.  The 
width  of  the  flow  field  is  unknown  because  it  depends  on  the  water  surface  elevation. 

Simultaneous  satisfaction  of  all  variables  in  the  steady  form  of  the  shallow 
water  equations  for  specific  geometry  and  flow  rate  over  a  particular  portion  of  a 
channel  requires  use  of  iterative  procedures  because  the  equations  are  nonlinear.  The 
success  of  an  iteration  technique  depends  on  the  quality  of  the  initial  estimate.  An 
adequate  initial  estimate  of  the  flow  variables  is  difficult  to  establish  for  flow  fields 
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which  contain  features  such  as  standing  waves  and  hydraulic  jumps.  The  numerical 
model  can  use  a  time  marching  procedure  to  address  transient  processes  or  to  calculate 
a  steady  state  flow  condition  as  a  time  asymptotic  limit  of  a  time  dependent  problem. 
Therefore,  although  steady  flow  solutions  are  addressed  in  this  study,  side  boundary 
locations  consistent  with  initial  conditions  which  are  not  the  steady  state  solution  are 
assumed.  Given  these  assumed  initial  conditions  and  side  boundary  locations,  solution 
of  the  transient  shallow  water  equations  is  proceeded  in  time  until  the  steady  state 
solution  is  obtained.  This  process  is  known  as  "time  stepping"  to  steady  state  and  can 
be  viewed  as  an  iteration  procedure  in  time.  At  each  time  step,  Newton-Raphson 
iterations  are  used  to  solve  the  nonlinear  equations.  As  the  computed  flow  field 
evolves  from  the  specified  assumed  initial  flow  conditions  and  initial  free 
surface/channel  boundary  location  to  the  steady  state,  the  side  boundaries  are  adjusted. 
This  constitutes  a  moving  boundary  problem. 

Objectives 

The  primary  focus  of  this  study  is  the  development  of  a  method  for  analyzing 
flow  in  trapezoidal  high-velocity  channels  where  the  flow  can  be  subcritical  or 
supercritical  and  may  undergo  transition  from  one  regime  to  the  other.  This  analysis 
is  primarily  concerned  with  water  surface  elevation  predictions  because  the 
determination  of  the  water  surface  characteristics  is  essential  in  the  successful  design 
and  evaluation  of  flood  control  channels.  Flow  simulations  are  conducted  using  a 
depth-averaged  two-dimensional  finite  element  model  with  moving  boundaries.  When 
applied  to  supercritical  flow,  the  numerical  magnitude  of  the  nonlinear  terms  in  the 
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governing  equations  are  large.  Therefore,  a  proper  moving  boundary  algorithm  for 
supercritical  flow  must  be  designed  to  circumvent  computational  instabilities.  To 
avoid  computational  instabilities  a  novel  means  for  solving  the  boundary  displacement 
and  the  flow  variables  simultaneously  for  the  moving  boundary  problem  is  developed. 
Tests  of  the  computational  model  consisting  of  comparisons  with  laboratory  flume  data 
demonstrate  that  supercritical  flow  in  channels  having  sloping  sidewalls  can  be 
modeled  satisfactorily  using  the  shallow  water  equations. 

Specifically,  this  dissertation  includes  the  following: 

(a)  Extension  of  the  shallow  water  equations  to  include  a  moving  reference  frame. 

These  depth-averaged  two-dimensional  flow  equations  describe  the  effects  of 
moving  boundaries. 

(b)  Development  of  a  numerical  model  of  these  moving  boundary  shallow  water 

equations.  This  numerical  model  emphasizes  numerical  stability  for  advection 
dominated  flow  containing  shocks  such  as  oblique  standing  waves, 
characteristic  of  supercritical  flow. 

(c)  Demonstration  that  the  model  reproduces  an  analytic  solution  of  the  shallow  water 

equations  of  steady  ideal  flow  in  a  bend  of  a  V-shaped  channel. 

(d)  Model  tests  comparing  simulation  results  with  measurements  obtained  in 

laboratory  flumes.  The  tests  include  supercritical  flow  in  four  geometric 
configurations: 

(i)  a  trapezoidal-to-rectangular  transition, 

(ii)  a  trapezoidal  channel  having  a  horizontal  curve. 
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(iii)  a  rectangular-to-trapezoidal  transition,  and 

(iv)  a  trapezoidal  channel  having  flow  obstructions  geometrically 

similar  to  bridge  piers. 

Outline  of  Dissertation 

A  description  of  flow  and  design  features  of  trapezoidal  high-velocity  channels 
has  been  given  above.  Chapter  2  reviews  the  methods  that  other  researchers  have 
proposed  for  modeling  moving  boundary  flow  problems.  In  Chapter  3  the  depth- 
averaged  two-dimensional  conservation  equations  are  developed  to  include  the  effects 
of  a  moving  reference  frame.  The  equations  are  established  in  this  form  to  enable  the 
formulation  of  a  moving  finite  element  model  for  advection  dominated  flow  fields 
(Chapter  4).  Chapter  5  compares  computed  numerical  model  results  with  an  analytic 
solution.  Additional  model  testing  is  described  in  Chapter  6.  Chapter  7  presents  the 
conclusions  and  an  outline  of  additional  research  for  model  extensions. 


CHAPTER  2 


PREVIOUS  RESEARCH 

This  chapter  presents  a  review  of  previous  research  in  numerical  modeling  of 
moving  boundary  problems.  This  review  is  divided  into  three  main  parts  beginning 
with  a  brief  listing  of  physical  problems  involving  moving  boundaries.  Next,  some 
examples  of  efforts  directed  toward  solving  laterally-averaged  two-dimensional 
problems  are  given.  The  final  portion  of  this  chapter  discusses  various  methods  used 
to  model  depth-averaged  two-dimensional  flow  problems  having  moving  boundaries. 
These  modeling  schemes  can  be  divided  into  three  categories:  the  wetting  and  drying 
method;  the  transformation  method;  and  the  method  selected  for  this  study,  the  moving 
finite  element  method. 

Moving  Boundary  Flow  Problems 

Physical  phenomena  involving  moving  boundary  flow  problems  include  the 
solid-liquid  interface  associated  with  the  melting  and  solidification  of  metals,  glass, 
and  semiconductor  crystals  (Yoo  and  Rubinsky  1986,  McLay  1988)  and  slide  coating 
in  the  manufacture  of  photographic  products  (Kistler  and  Scriven  1984,  Christodoulou 
and  Scriven  1989).  A  review  of  the  early  efforts  at  modeling  these  problems  is  given 
by  Lynch  and  Gray  (1980). 

In  the  hydraulic  arena,  the  water  surface  of  open  channel  flow  is  a  free 
boundary  and  if  the  sidewalls  are  sloping,  the  side  boundaries  are  located  where  the 
free  surface  intersects  the  sloping  sidewalls.  Transient  three-dimensional  analysis  of 
open  channel  flow  in  trapezoidal  channels  must  consider  boundary  movement  in  the 
horizontal  and  vertical  directions.  Unsteady  two-dimensional  modeling  of  these  flow 
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fields  reduces  the  moving  boundary  to  a  free-surface  calculation  in  laterally-averaged 
models  or  to  a  problem  of  finding  the  waterline  in  depth-averaged  models. 

Laterally-Averaged  Models 

Free-surface  models  have  been  developed  for  flow  over  spillways  (Dcegawa  and 
Washizu  1973,  Diersch  and  Martin  1978,  and  Bettess  and  Bettess  1983),  flow  under 
sluice  gates  (McCorquodale  and  Li  1971),  unconfined  ground  water  flow  (Neuman  and 
Witherspoon  1971,  Finn  and  Varoglu  1976,  and  France  1980),  sloshing  flow  and 
solitary  wave  propagation/runup  (Betts  and  Hall  1978,  Ramaswamy  and  Kawahara 
1987,  Okamoto  et  al.  1992,  and  Takizawa  et  al.  1992),  flow  at  a  hydraulic  jump 
(Chippada  et  al.  1994),  and  dam  break  simulation  (Wang  and  Wang  1994).  This  list  is 
by  no  means  complete,  but  these  examples  are  a  good  representation  of  the  research 
effort  that  has  been  directed  toward  these  areas. 

Bettess  and  Bettess  (1983)  summarize  that  finite  element  modeling  of  free- 
surface  flow  typically  uses  one  of  the  following  approaches: 

(a)  Fix  the  element  mesh  and  then  vary  the  element  properties  to  simulate  the  actual 

boundary.  This  method  has  been  used  to  simulate  the  saturated  and  unsaturated 
zones  of  seepage  through  porous  media  where  the  kinetic  energy  of  the  flow  is 
small  (Anderson  and  Woessner  1992). 

(b)  Invert  the  problem  such  that  the  grid  coordinates  are  the  dependent  variables  and 

then  assuming  ideal  fluid  flow,  use  the  stream  function  and  velocity  potential 
as  the  independent  variables  (Finn  and  Varoglu  1976).  This  method  is  not 
applicable  to  flow  subject  to  boundary  resistance. 
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(c)  Construct  a  finite  element  mesh  from  the  channel  bed  to  an  assumed  free  surface 
and  then  move  the  mesh  in  accordance  with  iterative  solutions  of  the  free 
surface,  until  convergence  is  achieved  (Ikegawa  and  Washizu  1973,  Betts  1978, 
and  Washizu  et  al.  1976). 

Each  of  these  approaches  can  be  applied  to  problems  where  the  location  of  the 
boundary  is  unknown  a  priori.  However,  modeling  supercritical  flow  while  including 
boundary  frictional  effects  excludes  the  use  of  methods  (a)  and  (b).  The  ideas 
associated  with  method  (c)  are  similar  to  many  of  those  used  to  solve  the  depth- 
averaged  equations  applied  to  subcritical  flow  fields.  These  iterative  methods  of  mesh 
adjustment  for  the  shallow  water  equations  are  described  in  the  next  section. 

Depth-Averaged  Models 

A  common  method  of  modeling  open  channel  flow  with  sloping  sidewalls  uses 
a  fixed  grid  with  an  imaginary  fixed  vertical  barrier  near  the  boundary,  as  shown  in 
figure  2.1,  rather  than  a  moving  boundary.  However,  this  does  not  produce  the  correct 
incident  wave  reflection  due  to  frictional  differences  (Lynch  and  Gray,  1978).  That  is, 
the  bed  drag  that  would  occur  in  the  absence  of  the  imaginary  fixed  barrier  is  not 
reproduced.  Frictional  losses  near  the  waterline  in  the  real  system  are  not  modeled 
using  an  imaginary  fixed  barrier. 

This  method  could  be  applied  to  high-velocity  channels;  however,  in  the 
supercritical  flow  regime,  the  steady  state  waterline  may  have  significant  spatial 
variation  along  the  channel.  This,  in  conjunction  with  the  possibility  that  the  depth  at 
the  imaginary  barrier  may  become  zero  as  the  solution  is  time  stepped  to  steady  state. 
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Figure  2.1.  Half  section  sketch  showing  imaginary  fixed  vertical  barrier. 

requires  that  the  imaginary  barrier  be  offset  far  from  the  expected  final  waterline. 

Also,  since  the  celerity  of  gravity  waves  is  a  function  of  flow  depth,  the  plan  view 
shapes  of  standing  waves  in  the  real  system  are  not  properly  modeled  using  a  fixed 
vertical  barrier. 

Depth-averaged  moving  boundary  flow  problems  have  been  solved  using  one  of 
three  methods.  The  three  methods  are  the  wetting  and  drying  method,  the 
transformation  method,  and  the  method  of  moving  finite  elements. 

Wetting  and  Drying  Method 

The  most  commonly  used  method  employs  a  fixed  spatial  grid  with  the  grid 
extended  to  cover  the  boundary  motion.  Subsequently,  the  grid  is  adjusted  iteratively 
to  cover  the  domain  width  as  a  function  of  the  computed  depth. 

This  method  typically  is  implemented  using  wet  and  dry  nodes.  If  the  nodes  in 
a  grid  cell  are  dry,  then  the  cell  is  "turned  off".  The  cell  is  no  longer  considered  in 
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the  flow  domain.  This  method  has  been  applied  to  the  two-dimensional  shallow  water 
equations  using  both  the  finite  difference  (Reid  and  Bodine  1968,  Leendertse  1970) 
and  finite  element  methods  (Holz  and  Withum  1977).  Lynch  and  Gray  (1980),  discuss 
the  early  attempts  at  modeling  moving  boundary  problems  in  this  manner.  More 
recent  finite  element  solutions  of  the  shallow  water  equations  using  wetting  and  drying 
are  given  by  Kawahara  and  Umetsu  (1986)  and  Leclerc  et  al.  (1990). 

Accurate  solutions  of  the  shallow  water  equations  are  difficult  to  obtain  using 
the  wetting  and  drying  method  because  this  method  does  not  conserve  mass  or 
momentum.  Each  time  a  cell  is  added  or  subtracted  to  the  domain  there  is  an 
instantaneous  change  in  the  momentum  near  the  boundary  and  in  the  total  volume  of 
fluid  in  the  flow  field.  Holz  and  Nitsche  (1980)  propose  a  scheme  to  insure  mass 
conservation;  however,  their  model  does  not  conserve  momentum  which  requires  the 
dynamics  on  all  partly  wetted  elements  be  considered.  The  primary  concerns  of  using 
the  wetting  and  drying  concept  to  model  flow  in  trapezoidal  high-velocity  channels  are 
that  the  method  can  produce  a  jagged  grid  rather  than  the  smooth  boundary  which 
actually  occurs,  rapid  changes  in  boundary  location  followed  by  periods  of  no 
boundary  movement,  and  stability  problems  in  advection  dominated  flow  (discussed  in 
Chapter  4).  In  supercritical  flow,  the  jagged  boundary  itself  will  create  shocks  and 
steady  state  convergence  may  not  be  obtainable  without  using  a  very  fine  grid  near  the 
boundaries. 

Transformation  Method 

A  second  method  of  simulation  consists  of  a  transformation  procedure  whereby 
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the  grid  is  transformed  back  in  time  to  its  initial  configuration.  Lan  et  al.  (1992) 
apply  this  method  to  the  one-dimensional  shallow  water  equations  to  compute  the 
water  surface  elevation  and  transverse  velocity  of  water  sloshing  from  side  to  side  in  a 
canal  having  a  parabolic  bathymetry  in  which  the  transverse  direction  is  aligned  with 
the  X  axis.  The  advantage  of  this  method  is  that  established  numerical  methods  for 
fixed-grid  problems  can  be  used  and  that  the  boundary  movement  has  no  limitation. 
However,  the  interpolation  of  the  bathymetry  relies  on  a  prescribed  spatial  relation. 

The  bathymetry  of  the  one-dimensional  example  given  by  Lan  et  al.  is  easily  described 
in  functional  form  because  the  bed  elevation  in  the  direction  of  flow  which  is  across 
the  channel  is  a  parabola  varying  only  in  the  x  direction.  Complex  boundary 
configurations  are  difficult  to  describe  in  functional  Euclidian  geometrical  form. 

The  flow  conditions  in  high-velocity  channels  are  determined  predominately  by 
the  geometric  configuration  of  the  channel  boundaries.  Therefore,  modeling  high- 
velocity  channels  requires  accurate  depictions  of  the  bottom  and  sidewall  geometry. 

The  transformation  scheme  seems  feasible  for  problems  with  simple  relations  between 
bottom  geometry  and  horizontal  position.  However,  it  does  not  seem  advantageous  for 
two-dimensional  modeling  of  high-velocity  channels  because  of  the  difficulty  in 
prescribing  the  channel  boundaries  in  a  three-dimensional  functional  form.  Three- 
dimensional  descriptions  of  complex  trapezoidal  high-velocity  channel  configurations 
functionally  must  describe  channel  side  slopes,  location  of  side  slope  toes,  and 
sidewall/side  slope  interception  in  wedge  type  rectangular-to-trapezoidal  (and 
trapezoidal- to-rectangular)  transitions.  Also,  plan  view  variations  such  as  horizontal 
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curves  and  transitions  are  difficult  to  describe  in  a  functional  form. 

Moving  Finite  Element  Method 

The  third  general  method  of  modeling  moving  boundary  problems,  using  finite 
elements,  involves  adapting  the  grid  to  adjust  the  size  and  number  of  elements  near 
the  side  of  the  domain,  coincident  with  the  solution  depth.  This  is  a  form  of  the 
moving  finite  element  method.  Papers  concerning  the  general  subject  of  moving  finite 
elements  have  been  written  by  Lynch  (1982)  with  application  to  phase  change 
problems  and  by  Miller  and  Miller  (1981)  and  Miller  (1981)  in  accurately  solving  the 
Burgers’  equation  (see  Finite  Element  Formulation  section  in  Chapter  4)  by 
concentrating  the  grid  density  near  flow  shocks.  A  thorough  description  of  the 
techniques  involved  in  implementing  the  method  of  moving  finite  elements  is  given  by 

Baines  (1994). 

There  are  two  common  means  of  accounting  for  flow  problems  in  which  the 
flow  varies  in  time.  Each  of  these  moving  grid  techniques  involves  a  transformation 
between  two  grids  whereby  nodal  points  are  moved  into  new  positions.  These  moving 
finite  element  methods  are  applicable  to  problems  that  are  time  stepped  from  assumed 
initial  conditions  to  the  steady  state  solution.  Both  of  these  methods  differ  from  their 
fixed  domain  counterparts  with  respect  to  temporal  derivatives.  One  method  describes 
temporal  derivatives  relative  to  a  fixed  coordinate  system.  Using  this  method,  the 
variation  in  time  of  the  finite  element  basis  (or  interpolation)  functions  must  be 
accounted  for.  This  manner  of  describing  temporal  derivatives  has  been  employed  for 
simulations  of  estuary-sized  problems  using  the  shallow  water  equations  by  Lynch  and 
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Gray  (1978)  and  using  the  shallow  water  wave  equations  by  Lynch  and  Gray  (1980) 
and  Siden  and  Lynch  (1988).  Lynch  and  O’Neil  (1981)  solve  the  two-phase  Stefan 
problem  and  Neuman  and  Witherspoon  (1971)  compute  the  free  surface  of  unconfined 
groundwater  flow  accounting  for  temporal  derivatives  in  this  manner. 

The  alternative  method  expresses  the  governing  equations  in  a  coordinate 
system  which  moves  with  respect  to  the  original  fixed  reference  frame.  Examples  of 
this  application  are  given  by  Ramaswamy  and  Kawahara  (1987)  using  the  arbitrary 
Lagrangian-Eulerian  (ALE)  description  of  the  fluid  domain  applied  to  a  width- 
averaged  model,  MacKinnon  and  Carey  (1993)  modeling  thermal  ablation  and 
consolidation  in  porous  media,  and  Akanbi  and  Katopodes  (1988)  simulating  the 
propagation  of  flood  waves  on  initially  dry  land. 

Lynch  (1982)  shows  that  these  two  methods  produce  identical  weighted 
residual  forms  of  the  partial  differential  equation  to  be  solved.  Following  the 
explanation  given  by  Lynch  (1982),  a  partial  differential  equation 


^+LU=/  (2.1) 

at 

is  examined  for  which  a  solution  of  U  =  U(x,t)  is  sought  for  a  deforming  grid.  Here,  t 
is  time  and  the  operator  L  contains  the  unknown  function  U  and  its  spatial  derivatives. 
The  weighted  residual  form  on  a  fixed  grid  is 
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where  the  integral  is  over  the  spatial  domain  Q,  U  =  U(x,t)  is  the  approximate 
numerical  solution  of  U,  and  Wj  are  the  weighting  functions.  The  numerical  solution 
of  U  (U)  is  obtained  via  interpolation  as 


0  =  E<i.jUj 


(2.3) 


where  (l)j  are  the  basis  functions  and  j  is  the  node  location. 
The  temporal  derivative  of  U  is 
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Lynch  shows  that 
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where  V®  is  termed  by  Lynch  as  the  "elemental  velocity".  Substitution  of  2.5  into  2.4 
yields 
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Then  the  weighted  residual  form  is 
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Lynch  points  out  that  the  term  V®  •  VU  "corrects  the  time  derivatives  to  account  for 
element  deformation".  Employment  of  this  technique  in  which  the  temporal 
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derivatives  are  described  relative  to  a  fixed  coordinate  system  produces  a  basis 
function  that  is  time  dependent  (equation  2.5). 

An  alternative  approach  using  coordinate  transformations  utilizes  spatial 
coordinates  %  which  move  relative  to  a  fixed  reference  frame  whose  spatial 
coordinates  are  X.  The  coordinate  transformation  can  be  stated  as 

X  =  X(x,t)  (2-8) 


and  the  time  derivative  of  a  dependent  variable  U  is  given  as 
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where  V  refers  to  the  fixed  reference  frame  coordinates,  X,  with  constant  t. 
Substitution  into  the  partial  differential  equation,  2.1,  yields 
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If  the  moving  coordinate  system,  %,  is  identified  with  the  local  coordinate 
system  which  is  also  used  for  the  finite  element  basis  functions,  then  the  numerical 
representations  of  the  temporal  derivatives  in  2.9  are 
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The  weighted  residual  form  of  2.10  is 
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which  is  identical  to  2.7,  yet  because  the  spatial  coordinate  system  which  moves 
relative  to  a  fixed  reference  frame  is  also  the  local  coordinate  system  used  to  define 
the  finite  element  basis  functions,  the  basis  functions  are  time  invariant. 

The  equality  of  2.7  and  2.13  illustrates  that  these  two  methods  of  describing 
moving  finite  elements  produce  the  same  weighted  residual  formulation.  One  can 
choose  to  express  temporal  derivatives  relative  to  a  fixed  coordinate  system  and  then 
account  for  the  variation  in  time  of  the  basis  functions  or  express  the  governing 
equations  in  a  coordinate  system  that  moves  with  respect  to  the  original  reference 
frame.  The  latter  method  is  used  in  this  dissertation.  The  reasoning  associated  with 
this  choice  of  methods  is  provided  in  the  Moving  Reference  Frame  Effects  section  of 
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CHAPTER  3 


EQUATIONS  FOR  MOVING  BOUNDARY  FLOW 

The  equations  used  to  model  open  channel  flow  are  presented  in  this  chapter. 
First  the  depth-averaged  two-dimensional  conservation  equations  are  given.  These 
equations  are  then  extended  to  include  the  effects  of  a  moving  reference  frame.  The 
equations  are  cast  in  this  form  to  describe  shallow  water  flow  whose  planform  domain 
varies  in  time. 

Shallow  Water  Equations 

Water  flow  in  open  channels  is  often  described  using  the  depth-averaged 
equations  of  mass  and  momentum  conservation,  commonly  referred  to  as  the  shallow 
water  equations.  The  shallow  water  equations  model  horizontal  motion  of 
incompressible,  Newtonian  fluid  flow  where  the  vertical  accelerations  are  negligible. 
The  dependent  variables  of  the  two-dimensional  fluid  motion  are  defined  by  the  flow 
depth  (h),  the  x-direction  component  of  depth-averaged  velocity  (u),  and  the  y- 
direction  component  of  depth-averaged  velocity  (v).  These  variables  are  functions  of 
the  independent  variables  x  and  y,  the  two  space  directions,  and  time  (t).  The 
variables  are  depicted  schematically  in  figure  3.1. 

Neglecting  free-surface  stresses  and  the  Coriolis  forces  as  these  are 
insignificant  in  high-velocity  channels,  the  shallow  water  equations  in  conservative 
form  are  given  as  (Liggett  1975,  Abbott  1979,  and  Praagman  1979) 

at  ax  ^  ^ 

for  the  conservation  of  mass.  The  conservation  of  momentum  in  the  x  direction  and  y 
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CHANNEL  BED 


Figure  3.1.  Flow  variables. 


direction  are  given  respectively  as 
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where 

Zq  =  channel  invert  (bed)  elevation 
g  =  acceleration  due  to  gravity 
n  =  Manning’s  roughness  coefficient 
p  =  fluid  density 

Cq  =  Cn/h^^^  is  a  dimensional  constant  (1  for  SI  units,  1.486  for 
English  units)  where  C  is  the  Chezy  coefficient 
^xx’  ^xy>  ^yx’  ^yy  Reynolds  stresses  per  unit  area  where  the  first 

subscript  indicates  the  direction  of  the  stress  and  the  second  indicates  the  axis 
direction  normal  to  the  face  upon  which  the  stress  acts.  The  depth-averaged  velocities 
u  and  V  are  temporally  averaged  values.  The  Reynolds  stresses  may  be  determined 
using  the  Boussinesq  approach  relating  stress  to  the  gradient  in  the  mean  currents 

=  2pv,-^  (3.4) 


and 
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ayy  =  2pv,|l  (3.6) 

where  Vj  is  the  turbulent  kinematic  viscosity.  The  turbulent  kinematic  viscosity  can  be 
approximated  empirically  as  a  function  of  depth  and  velocity  as  discussed  in  Chapter 
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The  depth-averaged  equations,  3.2  and  3.3,  are  derived  by  integration  of  the 
equations  of  motion  over  water  depth.  Depth  integration  of  the  momentum  equations 
results  in  momentum  transfer  terms  which  can  be  cast  in  forms  similar  to  stress  terms 
(Praagman  1979).  The  resulting  terms  are  commonly  referred  to  as  effective  stresses 
resulting  from  the  vertical  variation  of  horizontal  velocity  (Kuipers  and  Vreugdenhil 
1973). 

These  effective  stresses  resulting  from  the  vertical  variation  of  horizontal 
velocity  are  not  directly  included  in  3.2  and  3.3.  Exclusion  of  these  terms  implies  the 
assumption  that  vertical  variations  in  the  horizontal  velocities  are  negligible.  Some 
analysts  account  for  these  effective  stresses  resulting  from  depth  integration  using 
momentum  correction  factors.  The  momentum  correction  factor  applied  to  straight 
channels  can  vary  from  1.2  for  a  parabolic  vertical  distribution  of  horizontal  velocity 
to  1.05  for  many  turbulent  flows  (Liggett  1975).  However,  because  in  general  their 
values  are  not  known,  empirical  relations  are  often  used  to  quantify  the  effective 
stresses.  Coefficients  used  in  these  relations  are  based  on  flow  measurements 
(equation  4.5)  and  so  incorporate  both  the  effects  of  the  Reynolds  stresses  and  the 
effective  stresses  attributable  to  depth  integration. 

Related  velocity  and  depth  products  in  the  shallow  water  equations  are 
redefined  in  terms  of  unit  discharges 
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q  =  vh  (3-8) 

where  p  is  the  volumetric  discharge  per  unit  width  in  the  x  direction  and  q  is  the 
volumetric  discharge  per  unit  width  in  the  y  direction. 

The  variables  p  and  q  are  substituted  into  the  shallow  water  equations. 

Because  mass  is  conserved,  these  functions  are  smoother  than  u  and  v  throughout 
regions  in  the  flow  field  where  the  spatial  gradients  of  h,  u,  and  v  are  large  such  as  at 
an  oblique  standing  wave  or  at  a  hydraulic  jump.  Modeling  of  flows  containing  such 
gradients  (shocks)  using  the  shallow  water  equations  requites  that  the  depth-averaged 
equations  be  cast  in  this  form  (Abbott  1979).  The  governing  equations  are  rewritten  in 
terms  of  p  and  q  and  are  expressed  in  vector  form  as 
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Moving  Reference  Frame  Effects 

It  was  shown  in  Chapter  2  that  the  moving  finite  element  method  is  an 
appropriate  method  for  modeling  supercritical  flow  in  trapezoidal  channels.  The 
wetting  and  drying  method  does  not  conserve  mass  or  momentum  and  may  result  in 
unrealistic  jagged  waterlines  and  can  introduce  nonphysical  flow  shocks.  Although  the 
transformation  method  is  convenient  because  it  uses  existing  fixed  grid  technology,  it 
requires  functional  descriptions  of  the  channel  bed  which  are  difficult  or  impractical  to 
define  for  general  high-velocity  channels.  The  complex  boundaries  of  high-velocity 
channels  with  contractions,  expansions,  slope  breaks,  sidewall  slope  transitions,  curves, 
and  bridge  piers  would  be  very  difficult  to  describe  with  a  functional  geometrical 


32 


relation. 

It  was  shown  in  Chapter  2  that  two  moving  finite  element  methods,  one  using  a 
fixed  grid  with  a  temporally  varying  basis  function  and  the  other  expressing  the 
governing  equations  in  terms  of  a  coordinate  system  that  moves  relative  to  the  original 
fixed  reference  frame,  produce  the  same  weighted  residual  form  of  the  governing 
equations  (2.7  and  2.13).  The  latter  method  is  used  here  because  it  is  more  physically 
intuitive  and  modification  of  the  finite  element  basis  function  would  be  difficult  due  to 
the  particular  Petrov-Galerkin  test  function  employed.  The  governing  equations  are 
modified  to  account  for  the  moving  reference  frame  effects  rather  than  modifying  the 
finite  element  basis  function.  A  Petrov-Galerkin  test  function  that  is  composed  of  a 
Galerkin  part  and  a  non-Galerkin  component  (see  Petrov-Galerkin  Test  Function 
section  in  Chapter  4)  is  used  here.  Since  the  test  function  is  defined  in  terms  of  the 
basis  function,  modification  of  the  basis  function  to  account  for  a  moving  reference 
frame  is  likely  to  be  difficult  and  was  not  attempted. 

The  shallow  water  equations  describe  horizontal  motion  of  incompressible 
(Newtonian)  fluid  flow  relative  to  a  fixed  coordinate  system.  The  approach  used  here 
in  solving  the  moving  boundary  problem  transforms  the  governing  equations  derived 
for  a  fixed  coordinate  system  to  a  moving  coordinate  system  by  introducing  a 
coordinate  transformation  (MacKinnon  and  Carey  1993,  Akanbi  and  Katopodes  1988). 
The  shallow  water  equations  for  a  fixed  coordinate  system  are  expressed  in  terms  of  a 
moving  reference  frame  by  defining  a  new  set  of  independent  variables  (^.fij'C)  as  a 
moving  coordinate  system  which  is  related  to  the  fixed  reference  frame  (x,y,t)  such 
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that 


X  =  x(^,T|,X) 

(3.14) 

y  =  y(^,'n.'c) 

(3.15) 

and 


(3.16) 


The  time  derivative  of  the  flow  variables,  Q,  taken  in  the  moving  reference 
frame  is 


dQ  ^  dtdQ  ^  dx  dQ  ^  dy  dQ 
Bz  dz  dt  dz  dx  Bz  By 


where  BQ/Bt  is  the  time  derivative  in  the  fixed  coordinate  system  and  Bx/Bz  and  By/Bz 
are  the  moving  reference  frame  velocities  in  the  x  and  y  directions,  respectively.  The 
reference  frame  velocity  is  expressed  as 
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Using  the  chain  rule  of  differentiation. 
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The  shallow  water  equations  written  for  a  moving  reference  frame  are 
by  substitution  of  3.22  into  3.9 


dx 


dx 


l)  _  3(v,q)  ^ 

f  ^ 

3uj,  3vj. 

dy 

dx 


dy 


+  H  =  0 


Rearranging  yields 


^  +  A 

dx  dx 


X 


(  \ 


9uj.  3vj 


0  . 


(3.20) 


(3.21) 


(3.22) 


obtained 


(3.23) 


(3.24) 


+  H  + 
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The  vectors  F^*  and  Fy*  are  defined  as 


F/  =  F  -  u  O  = 


P  -  Ujh 


x  ■  ?■>" 


Aa., 


(3.25) 


and 


F  *  = 


-  F.  -  v,Q  = 


q  - 
pq 

^  -  v,p  -  -a,, 
■^  -  v,q  *  igh^  -  i. 


yy 


(3.26) 


The  shallow  water  equations  can  be  expressed  in  terms  of  a  moving  reference 
frame  as 


f 


dz  dx 


dF.*  ^Fy  „ 

+  _ L-  +  H  + 


ay 


IQ  .  0 


dx  dy 


(3.27) 


These  equations  can  be  used  to  describe  shallow  water  flow  whose  domain  may 


vary  in  time.  Provided  that  adequate  boundary  conditions  are  imposed  and  initial 
conditions  specified,  these  extended  shallow  water  equations  are  appropriate  for 
modeling  flow  in  channels  with  sloping  sidewalls. 


CHAPTER  4 


FINITE  ELEMENT  MODEL 

The  depth-averaged  two-dimensional  equations,  3.27,  expressed  in  terms  of  a 
moving  reference  frame  were  developed  in  the  previous  chapter.  These  equations, 
which  are  appropriate  for  modeling  depth-averaged  moving  boundary  flow  problems, 
form  a  system  of  mixed  nonlinear  elliptic  and  parabolic  partial  differential  equations. 
Although  these  equations  are  of  mixed  type,  in  applications  the  Reynolds  stresses  are 
relatively  small  and  the  system  behaves  much  as  hyperbolic  partial  differential 
equations.  The  Reynolds  stresses  are  important  however  in  modeling  flow  fields 
containing  eddies  associated  with  streamlines  separating  from  boundaries.  Such  flow 
fields  are  present  in  high-velocity  channels  in  the  vicinity  of  flow  obstructions  such  as 
bridge  piers. 

These  equations  in  general  cannot  be  solved  analytically.  Therefore  a 
numerical  solution  is  required.  The  finite  element  method  is  chosen  for  the  numerical 
representation  of  the  governing  equations  because  complex  geometric  details  of  high- 
velocity  channels  are  easily  represented  by  such  elements  and  boundary  conditions  are 
readily  imposed. 

This  chapter  begins  by  expressing  the  partial  differential  equations  in 
variational  form.  The  reference  frame  velocities  are  then  transformed  to  discrete  finite 
element  grid  velocities  (i.e.  nodal  velocities).  Then  the  finite  element  form  of  the 
governing  equations  is  presented.  Next,  to  simplify  boundary  condition  specifications 
and  reduce  functional  continuity  requirements,  a  weak  form  of  the  equations  is 
developed.  The  Petrov-Galerkin  test  function  which  provides  numerical  stability  is 
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then  given.  Following  this  the  finite  difference  equation  used  to  approximate  the 
temporal  derivatives  is  presented  and  then  the  appropriate  number  and  type  of 
boundary  conditions  are  examined.  The  movement  of  interior  nodes  used  to  reduce 
the  formation  of  ill  shaped  elements  during  boundary  nodal  displacement  is  then 
discussed  followed  by  methods  of  selecting  initial  conditions.  The  last  topic  of 
discussion  is  the  solution  procedure  which  is  presented  in  schematic  form.  Within  the 
solution  procedure,  the  Newton-Raphson  iterative  technique  is  used  to  solve  the 
nonlinear  finite  element  equations. 

Finite  Element  Formulation 

A  variational  formulation  of  the  governing  equations  involves  finding  a 
solution  of  the  dependent  variables,  Q,  using  the  test  function,  'F,  over  the  domain,  Q. 
The  variational  form  of  the  shallow  water  equations  related  to  a  moving  reference 
frame  is 


dx 


+  H  + 


faur 


UQ  =  0  . 


(4.1) 


The  finite  element  approach  used  is  a  Petrov-Galerkin  formulation  which 
incorporates  a  combination  of  the  Galerkin  test  function  and  a  non-Galerkin 
component  to  control  numerical  oscillations  due  to  advection.  Advection  dominated 
flow  refers  to  flow  conditions  in  which  the  Reynolds  number  defined  as 
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^  h\/u^  +  (4.2) 

V 

is  large.  Here,  v  is  the  molecular  kinematic  viscosity  of  the  fluid.  This  definition  of 
advection  dominated  flow  is  somewhat  ambiguous.  A  clearer  definition  from  a 
numerical  modeling  standpoint  is  that  advection  dominated  refers  to  flow  fields  which 
are  modeled  in  such  a  manner  that  the  mesh  Reynolds  number  (also  called  the  Peclet 
number)  is  large.  The  mesh  Reynolds  number,  defined  as 

Re  =  ^  (4.3) 

can  be  viewed  as  a  ratio  of  advection  to  diffusion  where  C  is  the  element  mesh 
parameter  (the  characteristic  length  of  the  grid  cells). 

Burgers’  equation,  which  contains  an  advection  and  a  diffusion  term,  serves  as 
an  analog  of  the  governing  equations.  The  linearized  Burgers’  equation  can  be 
expressed  as 


3u 


Un 


3u 


=  V, 


(4.4) 


at  "ax  ^3x2 

where  Uq  is  a  constant  advection  speed.  Stability  analysis  of  the  linearized  Burgers’ 
equation  shows  that  Galerkin  finite  element  schemes  using  linear  interpolation  (or 
finite  difference  methods  using  central  differences)  are  stable  for  mesh  Reynolds 
numbers  less  than  or  equal  to  2  (Anderson  ef  al.  1984).  Therefore,  in  the  work 
presented  here  the  term  advection  dominated  is  applied  to  modeled  flow  having  a 
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mesh  Reynolds  number  greater  than  2. 

An  order  of  magnitude  analysis  of  the  flows  simulated  and  numerical  grids 
used  here  provides  an  estimate  of  the  magnitude  of  the  mesh  Reynolds  number.  In 
open  channel  flow  where  the  turbulence  is  mainly  generated  by  bed  drag,  it  is 
reasonable  to  relate  the  horizontal  diffusion  of  momentum  to  the  bed  stress. 
Therefore,  the  turbulent  kinematic  viscosity  here  is  described  empirically  (Rodi  1980, 
and  Chapman  and  Kuo  1985)  as 

Vj  =  h  +  v^  . 

Here  f  is  the  Darcy-Weisbach  friction  factor  and  Cj,  is  a  coefficient  which  varies 
between  0.1  and  1.0.  It  has  been  shown  that  values  of  the  turbulent  kinematic 
viscosity  coefficient,  C^,  in  the  range  of  0.1  to  0.2  are  appropriate  for  modeling  high- 
velocity  channel  flows  (Berger  1993  and  Stockstill  and  Berger  1994).  In  design 
practice  an  engineer  should  conduct  simulations  with  various  values  of  C^,  and 
Manning’s  n  using  this  range  to  establish  an  envelope  of  possible  resulting  flow 
conditions.  This  turbulence  model  which  indicates  variation  of  turbulent  kinematic 
viscosity  with  local  depth  and  velocity,  is  however  used  in  the  present  study  to 
determine  a  constant  turbulent  kinematic  viscosity  for  the  entire  flow  field  based  on 
depth  and  velocity  for  uniform  flow.  This  simple  turbulence  model  is  deemed 
adequate  for  the  purposes  of  this  dissertation  to  model  high-velocity  channels  where 
the  primary  concern  is  obtaining  an  accurate  solution  of  the  flow  depth  whereby  the 
sidewall  heights  can  be  appropriately  designed  to  contain  flood  flows. 
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Substituting  4.5  into  4.3  results  in  an  expression  of  the  mesh  Reynolds  number 
in  terms  of  variables  for  which  an  order  of  magnitude  (or  range  of  order)  is  known. 

Re,  =  — !_  (4.6) 

c.yTh 


The  mesh  characteristic  length  and  the  flow  depth  are  of  the  same  order  of 
magnitude, 

1  =  0(10®) 

h 

where  0( )  is  the  "order  of.  The  value  of  Cj,  varies  as 

0(10"^)  <  Cb  <  0(10®)  . 

The  Darcy-Weisbach  friction  factor  for  the  hydraulically  smooth  boundaries  of 
concrete-lined  high-velocity  channels  is  0(10'^).  The  mesh  Reynolds  number  is  seen 
to  vary  from  O(IO^)  for  a  Cb  value  of  0(10®)  to  0(10^)  for  Cb  of  0(10-^).  For  the 
problems  addressed  here,  the  mesh  Reynolds  number  is  always  greater  than  2,  the 
model  must  use  a  numerical  scheme  which  is  stable  in  advection  dominated  flow. 

In  the  finite  element  form,  the  reference  velocity  is  taken  to  be  the  velocity  of 
the  grid,  =  Vg  =  (Ug,Vg).  Therefore,  F^*  =  -  QUg  and  Fy*  =  Fy  -  QVg.  Then  the 

finite  element  form  of  3.27  is 
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dy 


for  each  i  . 


(4.7) 


The  symbol indicates  the  discrete  value  of  the  quantity,  e  indicates  a  particular 
element,  and  i  indicates  a  particular  test  function.  The  grid  velocity  is  determined  in 
terms  of  the  nodal  displacement  rate.  The  displacement  vector  of  a  node  is 

S  =  |s|6  =  (|s|e. ,  |s|e,)  H-8) 

where  Isl  is  the  magnitude  of  the  displacement  vector  and  0  is  the  unit  vector  in  the 
direction  of  allowable  nodal  movement.  Specification  of  the  direction  of  allowable 
nodal  movement  is  arbitrary  at  this  point  in  the  model  development.  Because  the 
lateral  width  of  the  flow  field  is  sought,  specifying  the  direction  of  movement  as  a 
function  of  the  side  slopes  is  a  logical  choice.  Therefore,  boundary  nodes  located  on 
the  channel  side  slopes  are  restricted  to  moving  in  the  direction  of  maximum  channel 
side  slope  (see  figure  4.1).  That  is,  §  is  a  unit  vector  in  the  x-y  plane  directed  toward 
the  maximum  side  slope.  Given  that  s  is  the  nodal  displacement  vector,  then  the  grid 
velocity  is  Vg  =  ds/dt. 

The  moving  coordinate  system  is  identified  with  the  local  coordinate  system 
which  is  also  used  for  the  basis  functions.  The  geometry  and  flow  variables  are 
subsequently  interpolated  using  the  Lagrange  basis  functions 


42 


Figure  4.1.  Definition  sketch  of  the  direction  of  allowable  nodal  movement,  6. 


Q  =  y  (|)  Q 

Y  '  '  (4-9) 

and 

s  =  X)  (4.10) 

j 

where  ^  is  the  basis  function  and  j  is  the  node  location.  Bilinear  triangular  and 
quadrilateral  elements  are  used  with  nodes  at  the  element  comers.  Figure  4.2  shows 
the  two  bilinear  elements  used  in  terms  of  local  coordinates  ^  and  rj. 

The  test  function  used  (to  be  elaborated  on  later  in  this  chapter)  is 

\|/.  =  (1).I  +  (p.  (4.11) 

where  ([);  is  the  Galerkin  part  of  the  test  function  (identical  to  the  basis  function),  (pj  is 
the  non-Galerkin  part,  and  I  is  the  identity  matrix. 

To  facilitate  the  specification  of  boundary  conditions,  a  weak  form  (see  page 
44)  of  4.7  is  developed.  The  practical  weak  form  which  is  equivalent  to  that  given  in 
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Figure  4.2,  Transformation  from  global  coordinates  to  local  (bilinear  element) 
coordinates. 


4.7  is  found  using  an  integration  by  parts  procedure.  The  weak  form  includes 
boundary  integrals  that  are  actually  natural  boundary  conditions  for  the  shallow  water 
equations.  Using  Gauss’s  Theorem,  integration  by  parts  of  the  terms 


and 


are  given  as 
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,3(t),  i.  , 

dQ^  +  f  (1,.FX  dr 


(4.12) 


and 


dr 


(4.13) 


where  (n^.Hy)  =  n  is  the  unit  vector  outward  normal  to  the  boundary  r. 

Incorporating  the  results  of  the  integration  by  parts  procedure  into  4.7  yields 
the  weak  form  of  the  equations  which  is 


A 


f 

( 

Y 

3u. 

3v 

H  + 

—1^—L 

Q 

V 

)- 

da 


+ 


for  each  i 


(4.14) 


where  ~  is  omitted  for  clarity;  the  variables  are  understood  to  be  discrete  values. 

This  form  of  the  equations  is  named  weak  because  the  functional  continuity 
requirements  of  the  vectors  Fj^  and  Fy  in  the  terms  (J)  3Fj(  /dx  and  <j)  3Fy*/3y  have  been 
weakened  .  The  weak  form  only  requires  that  the  basis  functions  be  continuous 
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whereas  the  original  form  requires  a  basis  function  with  continuous  first  derivatives. 
Also,  the  weak  formulation  essentially  relaxes  the  boundary  condition  enforcement 
such  that  natural  boundary  conditions  are  enforced  in  an  average  sense.  Treatment  of 
natural  boundary  conditions  is  discussed  in  the  Boundary  Conditions  section  of  this 
chapter. 

Equation  4.14  can  be  simplified  first  by  realizing  that  =  F^*(Q,s)  and  Fy  = 
Fy*(Q,s)  and  expanding  the  momentum  matrix  terms  as 

=  (4.15) 

3x  0Q  dx  ds  dx 


and 


dFy  ^  dFy  aq  ^  (4.16) 

dy  dQ  dy  ds  dy 

Next,  the  convection  matrices  A*  and  B*  are  defined  as 

-Ug  1  0 

(c^-u^)  (2u-Ug)  0  (4-17) 

-(uv)  V  (u-Ug) 


. ,  <(Q)  . 
ag 


and 
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.  =  ^Fy‘(Q) 


-V,  0  1 

-(uv)  (v-Vg)  u 
(c^-v^)  0  (2v-vp 


(4.18) 


where  c  =  (gh)*'^  is  the  celerity  of  a  small  gravity  wave  in  shallow  water  resulting 
from  a  disturbance  in  the  local  depth.  The  convection  matrices  are  used  to  define  the 
characteristic  based  test  function.  A  discussion  including  the  development  of  the  test 
function  used  here  is  provided  in  the  next  section. 

Evaluation  of  the  derivatives  with  respect  to  s  gives 


and 


9s  9s 


p  -  u^h 


-u  D  +  -gh^  -  ha 


-  u  q  -  ha, 
h  * 


yx 


(4.19) 


9f; 

'dT 


9 

9s 


q  -  v^h 

f  -^sP  -  Ky 

-  V  q  +  — gh^  -  ha 

h  «  2^ 


yy 


9v 

Q 

9s 


(4.20) 


Therefore,  the  momentum  matrices  are 
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^'’x*  -  (4.21) 

ax  ^  ax  ax^ 

and 

^  -  ^O  (4.22) 

dy  ay  ay 

The  resulting  weak  form  which  is  found  by  substituting  4.21  and  4.22  into  4.14  is 
given  as 


au,  dv 

+  +  d) _ *  +  _ i 

^  ax  ay 


(4.23) 


♦  ^(fx  >  F;n,)dr  =  0. 

r 


The  final  simplification  is  made  by  introducing  a  notation  change.  A  new  vector  is 
defined  as 


(4.24) 


Reduction  of  4.23  produces  the  finite  element  description  of  the  shallow  water 
equations  modified  to  account  for  moving  reference  frame  effects 
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(4.25) 


If  the  grid  boundary  nodes  do  not  move,  that  is  if  Ug  =  Vg  =  0,  then  G  =  0,  F^* 
=  Fy*  =  Fy,  A*  =  dFJdQ,  and  B’  =  3Fy/3Q  and  the  model  reduces  to  a  standard 
fixed-grid  finite  element  formulation  as  it  should. 

Petrov-Galerkin  Test  Function 

The  particular  finite  element  test  function  employed  here  follows  closely  that 
first  given  by  Berger  (1992)  applied  to  the  shallow  water  equations  modified  to 
account  for  bed  curvature  effects  and  by  Berger  (1993)  and  Berger  and  Stockstill 
(1993)  for  application  to  advection  dominated  flows.  Each  of  these  studies  were 
directed  toward  fixed  grid  formulations  of  the  shallow  water  equations.  The  fixed  grid 
formulation  is  appropriate  for  modeling  bounded  flows  in  channels  having  vertical 
sidewalls.  Their  findings  show  that  the  test  function  which  is  weighted  upstream 
along  characteristics  provides  the  numerical  dissipation  needed  to  reduce  the  short 
wavelength  oscillations  associated  with  advection  dominated  flow.  It  is  well  known 
that  short  wavelength  oscillations  are  a  result  of  discretization.  Fourier  components  of 
the  equations  can  be  composed  of  all  wavelengths  as  a  result  of  boundary  conditions 
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and  computer  round-off  (Finder  and  Gray  1977).  Wavelengths  shorter  than  two  times 
the  grid  spacing  can  not  be  resolved  and  are  aliased  into  longer  wavelengths  (Abbott 
1979).  This  results  in  spurious  oscillations  if  no  numerical  dissipation  mechanism  is 
applied.  The  basic  idea  of  Berger’s  method  is  to  apply  diffusion  as  a  result  of  flow 
parameters.  The  scheme  is  a  form  of  the  Streamline  Upwind  Petrov-Galerkin  (SUPG) 
method  of  Hughes  and  Brooks  (1982).  Stockstill  and  Berger  (1994)  have  applied  this 
finite  element  model  to  various  rectangular  high-velocity  channel  configurations. 
However,  the  present  work  extends  Berger’s  test  function  to  account  for  moving  grid 
effects  so  that  flow  in  trapezoidal  high-velocity  channels  can  be  simulated. 

The  Petrov-Galerkin  part  of  the  test  function,  (pj,  is  defined  (Berger  1993,  and 
Berger  and  Stockstill  1993)  as 


9i  =  P 


(  ^ 
dd).  3(1); 

Ax — LA  +  Ay-_B 
dx  dy 


(4.26) 


where  P  is  a  dimensionless  number  between  0  and  0.5,  and  (])  is  the  linear  basis 
function.  This  part  of  the  test  function  serves  to  dissipate  spurious  oscillations  at  short 
wavelengths  attributed  to  advection  dominated  flows.  The  coefficient  |3  is  a 
dissipation  parameter.  Larger  values  of  P  result  in  large  dissipation  of  the  short 
wavelength  components  which  provides  stability,  yet  small  values  of  p  produce  more 
accurate  results  (Berger  1993).  In  4.26  Ax  and  Ay  are  representative  elemental 
dimensions  having  units  of  length.  Following  the  relation  given  by  Katopodes  (1984) 
the  elemental  lengths  are  chosen  as 
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Ax  =  2 


N 


dx 


^2  / 
+ 


(4.27) 


and 


Ay  =  2 


dy 


(4.28) 


where  |  and  T]  are  the  local  coordinates  defined  from  -1  to  1  (see  figure  4.2). 

A  and  6  are  computed  in  the  following  manner.  Beginning  with  the  definition 
of  the  convection  matrices  A*  and  B*  (4.17  and  4.18,  respectively),  the  characteristic 
matrices  can  be  rewritten  as 


A^  P  =  A* 


(4.29) 


and 


T  ‘  Ay  T  =  B- 


(4.30) 


where  A^  =  lA,  and  Ay  =  ly  are  the  matrices  of  the  eigenvalues  of  A*  and  B*, 
respectively.  The  matrices  P  and  P  ’  consist  of  the  right  and  left  eigenvectors  of  A* 
and  T  and  T  '*  are  made  up  of  the  right  and  left  eigenvectors  of  B*.  Then  A  and  ^ 
are  defined  as 


A  «  A 

A  s  p-i  A^  P 


(4.31) 


and 
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B  =  T  *  A.  T 


(4.32) 


where 


p-i  =  J_ 

2c 


0  1  -1 
0  (u+c)  -(u-c) 
1  V  -V 


(4.33) 


A.  = 


0 


0 


0 


*  7l 


0 


0 

0 

K 


+  7l 


(4.34) 


P  = 


-V  0  1 

-(u-c)  1  0 

-(u+c)  1  0 


(4.35) 


"p-i  —  ^ 

2c 


0  1 

1  u 


-1 

-u 


0  (v+c)  -(v-c) 


(4.36) 
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and 


(4.37) 


T 


-u  1  0 

-(v-c)  0  1 

-(v+c)  0  1 


(4.38) 


Here, 

=  u  -  Ug 

^2  =  (u  -  Ug)  +  c 

^3  =  (u  -  Ug)  -  C 

Yi  =  V  -Vg 

ll  =  (v  -Vg)  +  c 
Ys  =  (v  -Vg)  -  c 
and  c  =  (gh)^'^. 


The  test  function  employed  here,  V;,  has  two  parts,  a  Galerkin  portion  in  which 
the  test  function,  (|)j,  is  identical  to  the  basis  function  used,  and  a  Petrov-Galerkin 
portion,  (p^,  which  is  added  to  provide  the  numerical  stability  required  in  modeling 
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high-velocity  flows  containing  shocks  such  as  oblique  standing  waves  and  hydraulic 
jumps. 

Temporal  Derivatives 

The  semidiscrete  model  given  in  4.25  is  an  unsteady  form  of  the  governing 
equations.  Steady  state  solutions  are  obtained  by  marching  the  solution  in  time  until 
convergence  is  obtained.  Use  of  this  time  dependent  approach  requires  discretization 
of  the  temporal  derivatives  of  4.25.  Here  a  finite  difference  expression  is  used  for  the 
temporal  derivatives.  The  expression  for  the  temporal  derivative  of  the  variable 
vector,  Qj,  is 


( 

ag, 


=  a 


f  ^ 

Qm.l  _  Qm 


-  X® 


+  (1-a) 


_  Qm-1^ 


) 


X“  -  X”‘* 


(4.39) 


where  j  is  the  nodal  location  and  m  is  the  time  step.  An  a  equal  to  1  results  in  a  first 
order  backward  difference  approximation  and  an  a  equal  to  3/2  results  in  a  second 
order  backward  difference  approximation  of  the  temporal  derivative.  The  simulations 
conducted  here  are  all  for  steady  state  solutions,  therefore  an  a  =  1  is  used  exclusively 
because  time  accuracy  is  irrelevant.  This  eliminates  the  second  term  in  the  temporal 
derivative  (i.e.  only  two  time  levels  of  information  are  used)  thus  reducing  the 
computational  effort. 


Boundary  Conditions 

The  number  of  required  boundary  conditions  can  be  determined  by  examining 
the  characteristics  of  the  governing  equations.  Verboom  et  al.  (1982)  and  Drolet  and 
Gray  (1988)  show  from  the  theory  of  characteristics  that  the  number  of  boundary 
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conditions  required  is  equal  to  the  number  of  characteristic  half  planes  exterior  to  the 
domain  entering  the  domain.  Two  families  of  characteristics  exist  for  the  shallow 
water  equations.  The  first  family  is  a  cone  having  a  radius  equal  to  the  time-wave 
celerity  product,  t(gh)‘'^.  This  family  corresponds  to  the  concentric  rings  on  the  water 
surface  resulting  from  a  small  disturbance.  The  second  family  corresponding  to 
advection  is  the  plane  centered  on  the  axis  of  the  cone.  The  characteristics  are  paths 
on  which  information  about  the  flow  variables  travels.  Imposition  of  a  boundary 
prevents  information  from  traveling  in  the  same  manner  as  it  would  in  the  absence  of 
the  boundary.  Information  which  originates  within  the  domain  is  known  by  the  model 
and  therefore  requires  no  additional  effort.  However,  information  which  originates 
outside  the  domain  must  be  provided  to  the  model  as  boundary  conditions.  The 
required  number  of  boundary  conditions  for  the  inflow,  outflow,  and  sidewall 
boundaries  is  shown  in  table  4.1.  The  term  V„  is  the  outward  normal  component  of 
fluid  velocity. 

Open  Boundaries 

The  boundary  conditions  specified  at  inflow  boundary  nodes  for  subcritical 
inflow  are  u  and  v  or  p  and  q.  Supercritical  inflow  boundary  nodes  require 
specification  of  u,  v,  and  h  or  p,  q,  and  h.  The  tailwater  is  specified  as  the  only 
boundary  condition  for  subcritical  outflow  by  setting  the  depth  at  each  outflow 
boundary  node.  No  boundary  conditions  are  specified  for  supercritical  outflow. 

Side  Boundaries 

A  no  flux  boundary  condition  is  specified  for  sidewall  boundaries.  The  weak 
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form  of  the  governing  equations  includes  line  integrals  around  the  boundaries  (see  4.12 
and  4.13).  These  integrals  which  consist  of  flux  terms  are  the  system’s  natural 
boundary  conditions.  Natural  boundary  conditions  are  applied  to  all  moving 
boundaries  and  any  fixed  (vertical)  sidewall  boundaries.  These  boundaries  are  no  flux 
boundaries  in  that  is  there  is  no  net  flux  of  mass  or  momentum  through  them.  This 
boundary  condition  is  enforced  in  a  weighted  average  sense  through  the  weak 
statement.  Mass  flux  through  the  boundary  (  F  )  is  set  to  zero 

^  (p*  n^  +  q*  ny)dr  =  0  (4.40) 

r 


where  p*  =  p  -  Ugh  and  q*  =  q  -  Vgh.  Also,  there  is  no  net  momentum  flux  through  the 
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side  boundaries,  therefore 


f  [(up*)n^  +  (uq*)nj 


dr  =  0 


(4.41) 


and 


f  [(vp*)n^  +  (vq*)n 


dr  =  0 


(4.42) 


Sidewall  boundary  stresses  are  modeled  in  a  manner  agreeable  with  the  purpose 
of  the  work  developed  here.  The  purpose  is  to  provide  a  tool  for  evaluation  of  flood 
control  channels  in  which  determination  of  the  flow  depth  is  of  primary  importance. 
Because  applying  the  no  slip  condition  (i.e.  zero  velocity)  at  the  sidewalls  requires 
significantly  greater  grid  resolution  than  is  required  for  accurate  solution  of  the  flow 
depth,  sidewall  drag  is  treated  as  a  partial  slip  condition  by  approximating  the 
boundary  stress  terms  in  the  governing  equations  (integrated  along  the  sidewall)  by  the 
Manning  friction  loss  relationship 


/  dr  =  /  dr  (4.43, 


and 
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-  /  <l>i 


.ha  n)dr=  (  gql^LlZIZ.  . 
yy  yl  j  ^^^2  1,4/3 


(4.44) 


Therefore, 


f  <^^(FX  +  Fy-n^dr 

r 

r  0 


\jp^  + 

+  gp _ Ll _ 2_ 

n2  1,4/3 


C  ^ 

'^0 


]f] 

+  gq^- 


p2+  q2 


C  ^ 

'^0 


,4/3 


dr . 


(4.45) 


The  partial  slip  condition  allows  velocities  at  the  sidewalls.  If  the  depth  of 
flow  at  these  boundaries  is  set  equal  to  zero  as  it  exists  in  the  real  system,  there  are 
singularities  (i.e.  division  by  zero)  in  the  F^’,  Fy*,  and  H  vectors.  To  avoid  this 
problem,  depths  at  the  side  boundaries  are  set  to  a  small  constant  value,  ix  These  side 
boundary  depths  can  be  set  to  less  than  five  percent  of  the  centerline  depth  in  a 
trapezoidal  channel  as  demonstrated  in  Chapter  6  (see  figure  6.6).  This  introduces 
error,  but  if  the  specified  depth  at  the  boundary  is  sufficiently  small,  this  error  will  not 
significantly  affect  the  results. 

The  flow  velocity  and  grid  velocity  at  a  moving  boundary  are  related  by 

(V  -  Vg)  •  n  =  0  (4.46) 

where  n  is  the  unit  vector  outward  normal  to  the  moving  boundary.  The  boundary 
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conditions  at  the  moving  boundary  are 

=  h  (4.47) 

and 

X(x”)j  -  X(x“-i)j 

where  j  is  the  node  location  and  X  is  the  spatial  coordinate  vector.  The  time  step 
discretization  of  the  grid  velocity  is  identical  to  that  used  to  approximate  dQ/dx, 
Through  enforcement  of  the  natural  boundary  conditions  (equations  4.40,  4.41,  and 
4.42)  and  the  relationship  of  fluid  and  grid  velocities  (equation  4.46),  the  boundary 
locations  (waterlines)  are  calculated.  These  boundary  conditions  require  no  net  mass 

or  momentum  to  pass  in  or  out  of  these  side  boundaries  of  height  ^  at  the  waterline. 
The  no  flux  boundary  conditions  are  also  applied  to  the  vertical  sidewall  boundaries  of 
rectangular  channel  reaches.  Because  these  boundaries  do  not  move,  rectangular 
channel  reaches  are  modeled  using  a  fixed  grid  in  which  Vg  =  0. 

As  formulated  above,  there  are  potentially  four  degrees  of  freedom  at  each 
node:  hj,  pj,  qj,  and  Sj,  where  Sj  is  the  magnitude  of  displacement  of  node  j  in  the 

A 

direction  0  which  is  here  taken  as  the  direction  of  maximum  side  slope.  In  the  flow 
field  the  grid  does  not  move  and  so  p,  q,  and  h  are  the  unknowns.  However,  a  side 
boundary  on  a  sloping  sidewall  is  a  moving  boundary,  but  the  depth  of  flow  on  this 
boundary  is  constant  over  time  (h  =  ^  and  is  specified.  Rather  than  calculating  h  at 
this  moving  boundary,  the  boundary  location  and  the  boundary  velocity  must  be 


X(x^"i)j  -  X(x“)j 


+  (1-a) 


1 

I 
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calculated.  The  three  degrees  of  freedom  on  the  moving  boundaries  are  p,  q,  and  s. 
For  moving  boundary  nodes,  a  solution  is  sought  for  p^,  Qj,  and  Sj;  whereas,  the 
variables  to  be  solved  for  the  fixed  grid  nodes  are  pj,  qj,  and  hj. 

The  selection  of  boundary  flow  depth  values  is  problem  dependent.  The 
choice  of  an  appropriate  depth  at  the  moving  waterline  is  influenced  by  the  scale  of 
the  flow  field  simulated.  Choosing  ^  as  a  percentage  of  the  centerline  depth  given  for 
uniform  flow  is  a  logical  nondimensional  approach  for  trapezoidal  channels.  This 
method  is  used  in  the  simulations  reported  in  Chapters  5  and  6. 

Movement  of  Interior  Nodes 

The  boundary  conditions  associated  with  the  moving  boundary  nodes  are 
sufficient  for  describing  boundary  motion.  In  summary,  these  boundary  conditions  are 
that  there  is  no  net  flux  of  mass  or  momentum  through  the  moving  boundary  and  that 
the  depth  at  the  moving  boundary  is  a  constant.  However,  large  displacements  of  the 
moving  boundary  nodes  can  lead  to  element  shape  distortions.  Excessive  outward 
displacement  results  in  elements  adjacent  to  the  boundary  having  unacceptable  aspect 
ratios;  whereas,  inward  movement  can  lead  to  element  "folding"  (Baines  1994)  where 
nodes  move  across  element  lines  producing  ill  formed  elements  as  shown  in  figure  4.3. 

This  problem  can  be  overcome  by  manually  restructuring  the  computational 
mesh  every  few  time  steps  and  interpolating  from  the  results  of  the  last  time  step. 

This  is  quite  laborious  and  therefore  an  automatic  regridding  method  was  developed 
whereby  the  interior  nodes  lying  on  the  side  slopes  are  moved  each  time  step  as  a 
function  of  boundary  nodal  displacement.  Interior  side  slope  nodes  are  moved  along 


MOVING  NODE 


a.  Adjacent  elements  prior  to  nodal  displacement  (time  step  m) 


NODE  1  NODE  5 


NODE  3 

b.  Ill  formed  elements  resulting  from  excessive  nodal  displacement  (time  step  m+1) 
Figure  4.3.  Schematic  of  element  folding  from  excessive  nodal  displacement. 


grid  lines  (element  sides)  with  their  displacements  being  related  to  their  original 
relative  positions  across  the  side  slope.  For  example,  if  a  side  slope  is  evenly  divided 
into  4  elements  between  the  waterline  (where  the  exterior  node  is  a  moving  boundary 
node)  and  the  slope  toe  (where  the  interior  node  is  a  fixed  node),  then  given  the 
solution  of  the  boundary  node  displacement  for  a  time  step,  the  interior  side  slope 
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nodes  are  repositioned  at  the  quarter  points  across  the  side  slope.  MacKinnon  and 
Carey  (1993)  successfully  use  a  somewhat  similar  method  in  modeling  thermal 
ablation  and  consolidation  in  porous  media. 

After  new  x  and  y  coordinates  are  computed  for  the  moving  interior  nodes, 
flow  variables  and  the  invert  elevation  for  these  nodes  are  recomputed  using  a  linear 
interpolation  from  the  previously  computed  grid.  The  grid  velocities  of  these  moving 
interior  nodes  are  then  computed  as  a  displacement  rate  using  4.48  which  is  the  same 
equation  used  for  the  boundary  nodes. 

The  process  of  moving  interior  nodes  lying  on  the  side  slopes  as  a  function  of 
moving  boundary  nodal  displacement  helps  maintain  good  element  aspect  ratios  and 
reduces  element  folding.  By  only  allowing  the  nodes  associated  with  the  channel  side 
slopes  to  move,  the  invert  shape,  in  particular  the  side  slope  toe,  is  maintained. 

Initial  Conditions 

Initial  conditions  are  specified  at  the  beginning  of  each  simulation.  These 
initial  conditions  include  values  of  the  flow  variables  and  the  location  of  the  flow 
boundaries.  Because  these  initial  conditions  are  not  the  steady  state  solution  of  the 
two-dimensional  equations,  the  moving  finite  element  model  of  the  shallow  water 
equations  is  "time  stepped"  to  steady  state.  This  procedure  facilitates  the  calculation 
of  flow  fields  having  standing  waves  and  hydraulic  jumps. 

Initial  conditions  are  related  to  uniform  flow  conditions.  Given  the  discharge, 
cross  section  shape,  roughness,  and  slope  for  each  reach  of  the  channel,  the  Manning’s 
equation  is  used  to  compute  normal  depth  and  the  cross  sectionally  averaged  velocity 
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as  though  the  channel  were  prismatic.  The  depth  over  each  side  slope  is  then 
computed  from  the  cross  sectional  normal  depth.  The  velocity  computed  from  the 
Manning’s  equation  is  used  as  the  initial  velocity  for  straight  channels,  but  this  is 
cumbersome  for  geometrically  complex  channels  because  velocity  in  the  two- 
dimensional  model  is  a  vector  quantity.  Computing  the  x-  and  y-component  of 
velocity  for  channels  composed  of  numerous  reaches  having  various  cross  sectional 
and  plan  view  shapes  is  tedious.  Therefore,  the  initial  conditions  for  complicated 
channel  configurations  are  chosen  as  a  motionless  pool  at  normal  depth  for  each  reach. 
Initial  depths  and  velocities  are  specified  as  those  given  for  uniform  flow  for  straight 
channels;  whereas,  uniform  flow  depths  and  zero  velocity  are  given  as  initial 
conditions  for  geometrically  complex  channels. 

Solution  Procedure 

The  system  of  nonlinear  equations  is  solved  using  the  Newton-Raphson 
iterative  method.  Let  be  a  vector  of  the  nonlinear  equations  computed  using  an 
assumed  value  of  Lj  where 


f  '\ 

< 

Pi  “ 

Pj  ’ 

V 

Rj  is  the  residual  error  for  a 


if  node  j  is  an  interior  node 


if  node  j  is  a  moving  boundary  node  . 


(4.49) 


particular  test  function  i.  Subsequently,  Rj  is  forced 


toward  zero  as 
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17 


(4.50) 


where  k  is  the  iteration  number,  j  the  is  node  location,  and  dR^/dL.!"  is  the  Newton- 
Raphson  Jacobian  matrix.  This  linear  system  of  equations  is  solved  for  AL]^  and  then 
an  improved  estimate  for  L-'**  is  obtained  as 


=  L/  -  al;' 


(4.51) 


This  procedure  is  continued  until  convergence  to  an  acceptable  residual  error  is 
obtained. 

The  novel  method  for  solving  the  depth,  velocity,  and  domain  of  the  flow  field 
simultaneously  as  described  above  is  used  here.  A  simultaneous  solution  of  the  flow 
variables  and  the  boundary  displacement  circumvents  stability  and  convergence 


problems  arising  from  the  magnitudes  of  the  nonlinear  terms  in  the  governing 

equations.  The  partial  derivatives  of  the  residuals  with  respect  to  the  depth,  3RiV3hj, 

are  evaluated  for  the  fixed  grid  nodes,  whereas,  partial  derivatives  of  the  residuals  with 

respect  to  the  nodal  displacement,  3RiV3Sj,  are  evaluated  for  the  moving  nodes.  All 

derivatives  are  determined  analytically. 

A  contribution  of  this  study  is  the  determination  of  the  Newton-Raphson 

Jacobian  terms  involving  Sj  which  are  presented  in  Appendix  A. 

The  solution  proceeds  as  follows  for  each  time  step: 

(I)  Given  the  solution  from  the  last  time  step  or  initial  conditions,  two  time  levels  (m 
and  m+1)  of  information  is  at  hand  for  each  node  j: 
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hm+1  -m+l  ^  m+1  ^  m+1  _.m+l  \  m+1  /  \  m+1  urn  -m  m  j  m 

j  .  Pj  .  Qj  »  Xj  ,  Yj  ,  (Ug)j  ,  (Vg)j  ,  hj  ,  pj  ,  qj  ,  Xj  ,  ana  yj 

(note  for  initial  conditions,  the  variables  at  time  level  m  are  equal  to 
those  at  time  level  m+1) 


(II)  Time  Loop 


(A)  Update  variables: 

(1) 

u  m-l  _  U  HI 

hj  =hj 

then 

h” 

U  m+1 

(2) 

p”-‘  =  p” 

then 

Pj” 

=  Pj-"* 

(3) 

q,"*-!  =  q.™ 

then 

q™ 

= 

(4) 

Y™-!  -  V™ 

Xj  -  Xj 

then 

X."* 

li 

1 

(5) 

y”-l  =  y.™ 

then 

yj" 

=  yr 

(6) 

W  =  ( V 

m+1 

1 

(7) 

(vr  =  (vg) 

m+1 

j 

(B)  Newton-Raphson  Loop  (the  computed  variables  are  the  new  values  for 
time  level  m+1) 

(1)  Initialize  the  displacement  of  each  node  j: 

•  Sj  =  0 

(2)  Compute  Rj'^ 

(3)  Calculate  the  Jacobian  3Rf/3Lj‘': 

(a)  If  node  is  a  moving  boundary  node,  evaluate 


and 


aR,*^ 

■357' 


(b)  If  the  node  is  an  interior  node,  evaluate 

aR,*^  aR'=  aR*^ 

_ L,  — L,  and  -  . 

ah.  ap.  dq. 


(4)  Solve  the  linear  system  of  equations,  i.e.  compute  ALj'^. 

(a)  If  node  is  a  moving  boundary  node,  the  results  are 
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As, 


AL.  =  jAp. 

Kj 

(b)  If  the  node  is  an  interior  node,  the  results  are 


AL. 

j 


fAh.l 

} 

<  Ap.  > 


(5)  Update  the  variables  at  each  node  j: 

(a)  If  node  j  is  a  moving  node,  then: 


(i) 

Sj'^'^*  =  ASj^ 

(note  s/  is  initialized  to  zero  at 
the  beginning  of  each  iteration) 

(ii) 

p/">  =  p/  +  Ap/ 

(iii) 

q/"*  =  q/  +  Aq/ 

(iv) 

=  Xj*'  + 

(V) 

^  ^yj 

(Vi) 

2^*^  =  Zj‘‘  +  s/*‘(dzo/ds)j 

(vii) 

(Ug)/^*  =  (dx/dx)/^‘ 

(see  equation  4.48) 

(viii) 

(.\)r  =  (dy/dx)/"' 

(see  equation  4.48) 

where  (dVds)j  is  the  side  slope  in  direction 

of  §  =  (0^,  6y)  at  node  j. 

(b)  If  node  j  is  a  fixed  node,  then: 


(i)  h/^‘  =  h/  +  Ah/. 

(ii)  p/^i  =  p/  +  Ap/ 

(iii)  =  q/  +  Aq/ 

(6)  Evaluate  residual  error  using  the  max-norm  of  change  in  Froude 
number  from  the  last  iteration  (see  discussion  following  this 
solution  procedure  outline). 
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residual  error  =  max. 

1 


Frf'*  - 


Frf 


for  each  i . 


where 


+ 

11 

1 

(prf  -  (qff 

Frf  = 

(4  ^  (4 

44 

and  i  is  the  node  location. 

(a)  If  the  nondimensional  residual  error  is  greater  than  a 

specified  tolerance  (generally  of  the  order  10'^)  then 
return  to  step  (IIBl). 

(b)  If  the  nondimensional  residual  error  is  acceptable,  then 

move  the  interior  nodes 

(i)  move  interior  nodes  on  the  side  slope  along  an  element  side 

in  proportion  to  their  distance  from  the  associated 
moving  boundary  node 

(ii)  linearly  interpolate  the  bed  elevation  at  the  new  node 

position 

(iii)  linearly  interpolate  values  of  hj,  pj,  and  Qj  from  values 

computed  during  step  (IIB5) 

(iv)  compute  grid  velocities  of  moving  interior  nodes 

(Ug)j''‘^‘  =  (dx/dx)j‘'^‘  (see  equation  4.48) 

(Vg)j'‘^*  =  (dy/dT)j'''^‘  (see  equation  4.48) 

(III)  Advance  to  the  next  time  step,  go  to  step  (I). 
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The  Froude  number  change  is  a  nondimensional  representation  of  the  change  in 
the  values  of  the  flow  variables.  This  nondimensional  representation  avoids  problems 
of  scale  dependence.  The  same  convergence  criterion  is  applied  to  both  laboratory 
scale  problems  and  prototype  channels. 

The  coordinates  of  the  moving  nodes  are  updated  every  iteration  and  the 
Newton-Raphson  Jacobian  for  the  nodal  displacement  is  derived  accordingly 
(Appendix  A).  The  movement  of  interior  nodes  depends  on  the  displacement  of  the 
boundary  nodes.  However,  inclusion  of  these  effects  within  the  Newton-Raphson 
scheme  would  significantly  enlarge  the  coefficient  matrix  because  several  interior 
nodes  depend  on  each  moving  boundary  node.  To  avoid  this,  the  interior  nodes  are 
moved  once  each  time  step.  Therefore,  inclusion  of  the  effects  of  interior  node 
movement  within  the  Newton-Raphson  procedure  is  not  required. 


CHAPTER  5 


MODEL  VALIDATION 

This  chapter  provides  the  results  of  a  test  conducted  to  verify  that  the  model  is 
an  accurate  numerical  approximation  of  the  governing  equations.  It  is  important  to 
check  whether  the  discrete  numerical  model  is  actually  solving  the  shallow  water 
equations  properly.  This  verification  consists  of  comparing  model  results  to  an 
analytic  solution  of  the  shallow  water  equations.  The  particular  geometric 
configuration  chosen  for  this  test  is  a  horizontally  curved,  V-shaped  channel  of 
constant  radius  of  curvature  having  side  slopes  of  1  vertical  on  5  horizontal  and  no 
longitudinal  slope.  The  mild  side  slopes  were  chosen  to  accentuate  any  boundary  (i.e. 
waterline  location)  differences  between  the  model  and  the  analytical  results. 

Analytic  Equations 

An  analytical  solution  of  the  shallow  water  equations  written  in  cylindrical 
coordinates  is  determined  assuming  steady  irrotational  flow  of  an  inviscid  fluid  having 
no  radial  component  of  velocity.  The  z-direction  momentum  equation, 


U(0 


r  3ci) 


u. 


dz 


dP 


where 

r  =  the  radial  distance 

U(p  =  U(^(r)  =  the  co-direction  component  of  velocity 
u^  =  u^fr)  =  the  vertical  component  of  velocity 
p  =  the  fluid  density 
and 


(5.1) 
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P  =  the  pressure, 

under  the  hydrostatic  assumption  reduces  to 


ap 

dz 


=  -Pg  • 


Integrating  and  specifying  zero  pressure  at  the  surface, 


P(z)  =  pg(z  -  z) 


where  is  the  water-surface  elevation. 


The  r-direction  momentum  equation. 


A 


au,  au^  au, 

u _  + _ _  -  +  u _ - 

^  ar  r  act)  r  ^  dz  j 


ap 

"aF 


where  Uj.  is  the  radial  component  of  velocity,  reduces  to 


ap  ^U(o 

=  p - 

ar  r 


The  Bernoulli  equation  is  given  as 


-L.z.^=c, 

pg  2g 


where  Cj  is  a  constant.  Differentiation  of  5.6  with  respect  to  r  yields. 


(5.2) 

(5.3) 

(5.4) 

(5.5) 

(5.6) 

for  constant 


elevation  z 
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ap  _  ^(0 

17  ‘^“®"ar 


(5.7) 


Equating  5.5  and  5.7  and  integrating  provides  the  velocity  distribution  function  which 
is  identical  to  that  of  a  free  vortex. 


u 


(0 


(5.8) 


where  C2  is  a  constant. 

The  water-surface  gradient  found  by  substituting  5.3  into  5.5,  differentiating 
with  respect  to  r,  and  then  including  5.8  is  given  as 


^  2  ^2 

_  “(0  _  (5  9) 

gi^  ’ 

Test  Conditions 

The  normal  depth  concept  is  not  applicable  to  ideal  flow  in  a  channel  having 
no  longitudinal  slope.  Therefore,  initial  conditions  for  this  test  were  arbitrarily 
selected.  The  centerline  radius  and  depth  were  chosen  to  be  500.0  ft  and  5.0  ft, 
respectively.  Since  the  interest  of  this  research  is  on  modeling  supercritical  flow,  the 
lowest  local  Froude  number  chosen  was  1.0.  The  location  of  the  lowest  Froude 
number  is  at  the  channel  centerline.  Specifying  a  Froude  number  of  1.0  and  depth  of 
5.0  ft  results  in  a  centerline  velocity  of  12.69  ft/sec.  The  channel  side  slopes  were 
fixed  at  1  vertical  on  5  horizontal.  With  this  information,  the  velocity  distribution  and 
depth  variations  across  a  section  were  computed.  So  as  not  to  introduce  vorticity  into 
the  flow  field,  these  results  were  used  as  the  inflow  boundary  conditions  for  the  grid 
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nodes  located  at  O)  =  0  degrees.  The  initial  conditions  for  the  remainder  of  the  flow 
field  were  specified  as  a  uniform  velocity  distribution  and  a  level  water  surface. 

Inflow  boundary  and  initial  conditions  used  for  the  model  verification  are  shown  in 
figure  5.1.  Flow  at  the  outflow  boundary  was  supercritical  and  therefore  no  outflow 
boundary  condition  was  specified. 

The  model  parameters  used  are  presented  in  table  5.1.  The  boundary  depth,  ^ , 
was  chosen  to  be  5  percent  of  the  centerline  depth  resulting  in  steady  state  local 
Froude  numbers  of  4.7  at  the  inside  boundary  and  4.3  at  the  outside  boundary.  This 
choice  of  1i  resulted  in  an  error  of  approximately  0.25  percent  in  cross  sectional  area. 
This  error  was  deemed  acceptable  for  model  verification. 

The  curved  channel  was  discretized  into  6  elements  in  the  radial  direction  and 
the  length  of  the  channel  was  divided  into  207  elements  such  that  the  length-to-width 
ratio  of  each  element  was  approximately  2.  The  numerical  model  computational  mesh, 
shown  in  figure  5.2,  had  1456  nodes  and  1242  elements. 

Initial  time  step  selection  was  based  on  the  Courant  number  defined  as 


n.  -  +  \/gir 

{  . 
min 

“aT 


(5.10) 


where  is  the  minimum  grid  spacing  and  At  is  the  time  step.  The  numerical 
stability  constraint  for  many  explicit  models  of  hyperbolic  partial  differential  equations 
is  the  Courant-Friedrichs-Lewy  (CFL)  condition  (Anderson  et  al.  1984).  The  CFL 


a.  Inflow  boundary  conditions  (nodes  at  co  =  0  degrees) 


b.  Initial  conditions  (nodes  at  0  degrees  <  0)<  360  degrees) 

Figure  5.1.  V-shaped  channel  boundary  and  initial  conditions  (exaggerated  abscissa). 
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Table  5.1:  V-shaped  channel,  model  parameters. 


Condition 

Value 

a 

1.0 

p 

0.1 

n 

0 

V 

0 

a 

0.25  ft 

condition  limits  the  time  step  size  by  restricting  the  Courant  number  to  be  less  than  or 
equal  to  1.0.  The  model  developed  for  this  dissertation  is  fully  implicit  (equations 
4.39  and  4.48)  and  is  not  subject  to  the  CFL  condition;  however,  experience  in 
modeling  rectangular  high-velocity  channels  (Stockstill  and  Berger  1994)  has  shown 
that  a  time  step  corresponding  to  a  Courant  number  near  unity  is  a  good  choice  for 
startup. 

Simulations  were  conducted  using  a  time  step  of  1.0  sec  which  resulted  in 
Courant  numbers  of  1.1,  1.7,  and  1.0  at  the  inside  boundary,  channel  centerline,  and 
outside  boundary,  respectively.  The  model  converged  to  steady  state  in  250  sec. 

Test  Results 

The  analytic  solution  is  compared  to  model  results  at  the  cross  section  located 
at  an  CO  value  of  180  degrees.  The  flow  depths  are  compared  with  the  analytic 
solution  in  figure  5.3.  A  comparison  of  the  longitudinal  velocity,  u^^,  is  shown  in 
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figure  5.4.  The  largest  velocity  error,  (V^^dei  -  VanaiyticaiVVanaiyticai’  was  0.148  percent 
at  the  first  node  outside  the  centerline  node  where  the  simulation  velocity  at  this  node 
was  too  high. 


! 


Figure  5.2.  Numerical  model  computational  mesh  for  V-shaped  channel. 
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Figure  5.3.  Analytical  and  model  computed  flow  depths  for  V-shaped  channel 
(exaggerated  abscissa). 


RADIAL  DISTANCE  (r).  ft 


Figure  5.4.  Analytical  and  model  computed  longitudinal  velocity  for  V-shaped 


channel. 
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Since  there  are  no  energy  losses  throughout  the  flow  field,  the  specific  energy 
defined  as 


E  =  -^  +  h  +  z. 
2g  * 


(5.11) 


is  a  constant.  Choosing  the  bed  elevation,  Zq,  as  zero  at  the  channel  centerline,  values 
of  specific  energy  are  presented  in  table  5.2.  The  error  shown  in  table  5.2  is  the  ratio 
of  the  model  results  minus  the  analytic  solution  divided  by  the  analytic  solution.  The 
largest  error  in  the  computed  solution  was  0.051  percent  and  is  associated  with  the 
first  node  outside  the  centerline  node.  The  volumetric  flow  rate  computed  by  the 
model  through  the  cross  section  at  an  ®  value  of  180  degrees  differed  from  that 
specified  at  the  inflow  boundary  by  0.039  percent.  This  continuity  check  demonstrates 
that  the  model  conserves  mass.  These  results  show  that  the  model  is  a  numerical 


Table  5.2:  V-shaped  channel,  analytical  solution  and  model  computed  specific  energy. 


Radial  Distance,  ft 

Specific  Energy,  ft 

Percent  Error 

Analytic  Solution 

Model  Results 

477.4667 

7.5000 

7.4980 

-0.027 

484.9626 

7.5000 

7.4990 

-0.013 

492.5042 

7.5000 

7.5027 

0.036 

500.0000 

7.5000 

7.4993 

-0.009 

508.2810 

7.5000 

7.5038 

0.051 

516.6127 

7.5000 

7.4989 

-0.015 

524.8941 

7.5000 

7.5030 

0.040 
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representation  of  the  ideal  shallow  water  equations  with  moving  boundaries.  It  is 
noted  that  these  ideal  equations  exclude  boundary  friction  and  turbulence  effects. 


CHAPTER  6 


MODEL  TESTING 

Additional  model  evaluations  are  made  in  this  chapter  by  comparing  model 
results  with  laboratory  data.  These  evaluations  are  referred  to  as  model  testing  in  that 
tests  were  conducted  to  determine  the  shallow  water  model’s  ability  to  simulate  real 
flows. 

A  series  of  flume  tests  were  conducted  at  the  U.S.  Army  Engineer  Waterways 
Experiment  Station,  Hydraulics  Laboratory.  Four  geometric  configurations  were 
tested.  One  flume  consisted  of  a  trapezoidal  channel  having  a  horizontal  curve.  Three 
flume  tests  were  conducted  using  a  trapezoidal  channel  constructed  on  a  tilting 
platform.  These  tests  documented  the  flow  conditions  at  a  trapezoidal-to-rectangular 
transition,  a  rectangular-to-trapezoidal  transition,  and  flow  obstructions  geometrically 
similar  to  bridge  piers. 

The  flow  in  each  test  was  on  average  in  the  supercritical  regime.  The  cross 
sectional  average  Froude  number  given  as 

Fr  =  -pLr  (6.1) 


was  greater  than  unity.  Here  V  is  the  cross-sectional  averaged  velocity  (discharge 
divided  by  cross-sectional  flow  area  normal  to  the  flow  direction),  and  h^^^  is  the 
effective,  or  hydraulic  depth  (cross-sectional  flow  area  normal  to  the  flow  direction 
divided  by  the  free-surface  width).  However,  local  Froude  numbers  expressed  in 


terms  of  local  flow  variables  as 
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F  ^  l/u"  *  v"  (6,2) 

may  be  less  than  or  greater  than  unity  along  the  sloping  sidewalls.  The  flow  over  the 
sloping  sidewalls  may  be  subcritical  even  though  the  cross  section  on  average 
experiences  supercritical  flow. 

Frequently  the  water-surface  elevation  in  supercritical  channels  oscillates  in 
time  even  under  steady  inflow  boundary  conditions.  This  is  particularly  prevalent  in 
the  vicinity  of  shocks  such  as  oblique  standing  waves  and  hydraulic  jumps.  Water 
surface  oscillations  were  observed  in  each  of  the  flume  tests  conducted  in  this  study. 
The  flow  in  each  test  was  only  quasisteady  in  that  flow  depths  and  velocities  varied 
about  mean  values.  The  data  presented  herein  are  approximate  temporally  averaged 
values. 

Model  testing  described  in  this  chapter  was  conducted  in  a  progression  of 
difficulty  from  a  moving  finite  element  modeling  standpoint.  The  first  comparison 
was  with  supercritical  flow  in  a  trapezoidal-to-rectangular  transition.  Although  the 
first  test  was  not  a  demanding  test  for  the  moving  boundary  scheme,  it  did  show  the 
model’s  ability  to  simulate  oblique  standing  waves.  The  second  test  was  a  more 
general  two-dimensional  case  of  supercritical  flow  in  a  trapezoidal  channel  with  a 
horizontal  bend.  The  third  test  was  of  supercritical  flow  in  a  nectangular-to-trapezoidal 
transition  and  the  fourth  test  was  of  subcritical  and  supercritical  flow  in  a  trapezoidal 
channel  having  multiple  flow  obstructions.  The  flow  conditions  in  these  last  two  tests 
included  oblique  standing  waves  in  the  trapezoidal  channels. 
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The  model  parameters  used  in  each  of  these  tests  were  based  solely  on 
information  known  prior  to  model  simulations.  These  parameters  were  never  used  as 
model  calibration  devices  in  that  the  model  parameters  were  not  adjusted  to  improve 
the  simulation  results’  comparison  with  the  observed  laboratory  data. 

Trapezoidal-to-Rectangular  Transition 
Flume  tests  were  conducted  using  a  tilting  open  channel  flume  80  ft  long  and  1 
ft  high  with  a  base  width  of  2  ft  and  side  slopes  of  1  vertical  on  2.25  horizontal 
(figure  6.1).  The  flume  was  constructed  of  plastic  coated  plywood  covered  with  rolled 
epoxy  paint.  The  range  of  tilt  allowed  the  flume  to  be  set  on  hydraulically  steep 
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Figure  6.1.  Plan  and  elevation  views  of  tilting  flume. 
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slopes  which  produced  supercritical  flow.  Water  was  supplied  to  the  flume  by  a 
circulating  system.  The  water  was  pumped  from  a  sump  to  a  12-inch  supply  pipe. 

The  flume  discharge  was  returned  to  the  sump  through  a  drain  line.  The  discharge 
was  measured  using  a  factory  calibrated  12-inch  by  7.25-inch  venturi  meter  installed  in 
the  supply  line.  This  inflow  was  discharged  into  a  headbox  and  baffled  with  a  screen. 
Flow  at  the  downstream  end  of  the  channel  was  supercritical.  Boundary  conditions 
should  not  be  imposed  at  supercritical  outflow  boundaries  (table  4.1);  therefore,  no 
tailgate  was  used. 

Documentation  of  the  flow  conditions  consisted  of  photographs  and  measured 
flow  depths  and  velocities.  Flow  depths  were  measured  using  point  gages  with  a  limit 
of  reading  of  0.001  ft.  Velocity  measurements  were  made  with  a  commercial 
propeller  type  flow  meter  mounted  to  permit  measurement  of  current  from  any 
direction;  the  meter  elevation  was  also  adjustable.  The  0.46-inch  (outside  diameter) 
measuring  head  of  the  velocity  meter  was  constructed  of  a  five  bladed  PVC  rotor 
mounted  on  a  stainless  steel  spindle.  For  the  magnitude  of  velocities  reported  in  this 
study,  the  manufacturer  specification  lists  the  probe  accuracy  as  ±  0.05  ft/sec. 

A  trapezoidal-to-rectangular  transition  was  built  by  placing  vertical  walls  at  the 
top  of  the  side  slopes  starting  36.5  ft  from  the  upstream  end  of  the  flume  and  ending 
at  the  side  slope  toe  at  sta  50.  The  transition  rate  was  1  lateral  on  6  longitudinal.  The 
flume  was  rectangular  in  cross  section  from  sta  50  to  the  flume  end  (sta  80).  A  plan 
view  of  the  trapezoidal-to-rectangular  channel  transition  is  shown  in  figure  6.2.  Figure 
6.3  is  a  dry  bed  photograph  of  this  transition. 
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Figure  6.2.  Plan  view  and  sections  of  trapezoidal-to-rectangular  transition. 


Test  Conditions 

Model  testing  consisted  of  comparing  flume  data  obtained  with  a  discharge  of 
4.00  ft^/sec  on  a  slope  of  0.010  with  simulation  results.  The  test  setup  produced 
supercritical  flow  with  an  average  Froude  number  approaching  the  transition  (within 
the  trapezoidal  channel)  of  2.0.  This  initial  configuration  was  a  simple  test  for  the 
moving  boundary  model  because  the  straight  prismatic  trapezoidal  channel  flow 
conditions  were  nearly  uniform  with  smooth  waterlines  and  flow  boundaries  were 
established  by  the  vertical  sidewalls,  but  only  after  the  waterline  intersected  the 
vertical  transition  walls.  However,  as  the  supercritical  flow  over  the  sloping  sidewalls 
intercepted  the  vertical  transition  walls,  oblique  standing  waves  were  produced  within 
the  rectangular  portion  of  the  flume  as  shown  in  figure  6.4.  These  oblique  positive 
waves  were  generated  on  each  side  of  the  channel  and  met  at  the  channel  centerline 
where  the  maximum  local  flow  depth  was  produced.  The  standing  waves  were 
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Figure  6.4.  Flow  conditions  in  trapezoidal-to-rectangular  transition,  approach  Froude 
number  =  2.0  (looking  downstream). 
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reflected  from  the  sidewalls  for  quite  some  distance  downstream  of  the  transition. 
Although  this  first  case  was  not  a  definitive  test  for  the  moving  boundary  scheme,  it 
did  show  the  model’s  ability  to  simulate  oblique  standing  waves. 

Initial  mesh  construction  within  the  trapezoidal  reach  of  the  channel  required 
an  estimate  of  the  top  width  of  the  flow  so  the  location  of  the  side  boundaries  could 
be  established.  In  trapezoidal  channels  the  top  width  is  set  by  the  flow  depth, 
therefore  an  estimate  of  the  flow  depth  was  sought.  A  Manning’s  n  value  of  0.009  for 
the  smooth  boundaries  of  epoxy  painted  plywood  was  used  to  compute  normal  depth 
for  the  trapezoidal  channel  with  a  discharge  of  4.00  ft^/sec,  a  channel  slope  of  0.010,  a 
base  width  of  2.0  ft,  and  1  vertical  on  2.25  horizontal  side  slopes.  The  top  width  of 
the  initial  grid  was  determined  using  the  computed  normal  depth  applied  to  this 
channel  while  accounting  for  the  specified  depth  at  the  moving  boundary,  ti.  The 
lateral  horizontal  distance  from  the  side  slope  toe  to  the  waterline  was  set  equal  to  the 
difference  in  normal  depth  and  the  depth  at  the  waterline  divided  by  the  side  wall 
slope  (1/2.25).  The  initial  mesh  had  458  nodes  and  389  elements  and  is  shown  in 
figure  6.5.  The  initial  grid  contained  3  elements  across  the  channel  bottom  at  the 
upper  end  (headbox  outlet,  sta  2.82)  and  4  elements  across  the  bottom  from  upstream 
of  the  transition  (sta  28.8)  to  the  flume  end  (sta  80).  Side  slope  resolution  varied  from 
one  element  at  the  upper  end  to  2  elements  on  each  side  slope  in  the  vicinity  of  the 
transition  beginning  at  sta  30.0.  The  upstream  end  of  the  mesh  was  composed  of 
points  used  to  simplify  the  boundary  condition  input.  The  nodes  at  the  upstream  end 
of  the  channel  represented  a  2.82  ft-long  "headbox"  in  that  the  inflow  cross  section  (at 
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sta  0.0)  was  of  rectangular  shape  and  the  mesh  was  fixed  at  these  boundaries.  This 
"box"  transitioned  from  rectangular  to  the  actual  flume  shape  over  the  box  length  of 
2.82  ft.  Some  44  ft  upstream  of  the  transition  (waterline/vertical  wall  intersection) 
were  modeled  to  minimize  any  entrance  effects  caused  by  the  "headbox".  A  30-ft 
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Figure  6.5.  Numerical  model  computational  meshes  for  trapezoidal-to-rectangular 


transition. 
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long  reach  of  the  channel  downstream  of  the  transition  was  also  included  in  the 
calculations. 

The  inflow  boundary  was  supercritical  and  required  the  specification  of  p,  q, 
and  h  at  each  inflow  node.  Since  the  computational  mesh’s  x  axis  was  made  to  lie  on 
the  flume  centerline,  q  was  set  to  zero  and  each  value  of  p  was  determined  as  the  total 
discharge  divided  by  the  width  of  the  "headbox"  (3.215  ft).  Flow  at  the  outflow 
boundary  was  supercritical  in  the  rectangular  channel  and  so  no  boundary  conditions 
were  specified.  The  boundary  conditions  used  are  presented  in  table  6.1  and  the 
model  parameters  used  are  listed  in  table  6.2.  The  turbulent  kinematic  viscosity  was 
computed  (equation  4.5)  using  normal  depth  and  velocity  for  the  rectangular  reach  of 


Table  6.1:  Trapezoidal-to-rectangular  transition,  boundary  conditions. 


Boundary 

Flow  Regime 

Boundary  Condition 

Value 

Inflow 

Supercritical 

h 

0.280  ft 

P 

1.244  ft^/sec 

q 

0.0 

Outflow 

Supercritical 

h 

Not  Specified 

P 

Not  Specified 

q 

Not  Specified 
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Table  6.2:  Trapezoidal-to-rectangular  transition,  model  parameters. 


Condition 

Value 

a 

1.0 

p 

0.1 

n 

0.009 

0.1,  0.023  ft^/sec 

1^ 

0.01  ft 

channel.  This  was  the  major  area  of  interest  concerning  model  evaluation  because  this 
portion  of  the  channel  contained  the  oblique  standing  waves  generated  at  the  transition. 
The  chosen  parameters  were  never  adjusted  to  make  the  model  results  more  closely 
match  the  laboratory  data.  The  numerical  model  was  never  "calibrated"  to  replicate 
the  observed  quantities.  The  selected  constant  depth  (li  =  0.01  ft)  at  the  moving 
sidewall  boundaries  (figure  6.6),  was  less  than  5  percent  (3.6  percent)  of  the  depth  at 
the  channel  center  and  produced  a  cross  sectional  area  error  of  approximately  0.03 
percent  at  the  inflow  boundary.  Sidewall  boundaries  within  the  trapezoidal  portion  of 
the  channel  were  allowed  to  move  in  conjunction  with  the  depth  solution.  All  moving 
nodes  were  confined  to  move  in  the  direction  of  maximum  side  slope  with  the 
exception  of  each  sidewall  boundary  node  at  the  vertical  wall  where  the  trapezoidal-to- 
rectangular  transition  begins.  These  nodes  were  only  allowed  to  move  in  the  direction 
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parallel  to  the  transition  walls.  This  allowed  the  waterline  to  move  upstream  or 
downstream  as  the  solution  depth  increased  and  decreased,  respectively.  The 
remaining  sidewall  boundary  nodes  within  the  transition  and  all  nodes  within  the 
rectangular  portion  of  the  channel  were  fixed. 

Initial  conditions  for  these  simulations  were  the  flow  variables  given  for 
uniform  flow  in  each  of  the  two  typical  sections.  That  is,  normal  depth  (0.28  ft)  and 
velocity  (5.41  ft/sec)  for  the  trapezoidal  reach  were  used  as  the  initial  conditions  for 
the  upstream  reach  (sta  0.0  to  sta  46.36)  and  normal  depth  (0.34  ft)  and  velocity  (5.95 
ft/sec)  for  the  rectangular  reach  were  used  as  the  initial  conditions  for  the  lower  end  of 
the  flume  (sta  50  to  sta  80).  Approximate  mean  values  of  the  uniform  conditions  in 
the  upper  and  lower  reaches  were  used  for  the  initial  conditions  within  the  transition. 
Given  the  initial  and  boundary  conditions,  simulations  were  conducted  with  the  initial 
mesh  (figure  6.5a)  using  time  steps  of  0.1  sec  and  0.2  sec  for  a  total  time  of  13  sec. 

A  time  step  of  0.1  sec  resulted  in  a  Courant  number  of  approximately  0.9  at  the 
smallest  grid  element.  Although  the  model  is  fully  implicit,  the  computational  runs 
described  in  this  dissertation  typically  were  started  using  time  steps  corresponding  to 
Courant  numbers  near  1.0.  After  model  startup  the  time  steps  were  increased  to  speed 
up  steady  state  convergence. 

The  computational  mesh  was  refined  at  this  point  and  the  current  values  for  the 
coarse  mesh  were  interpolated  to  the  refined  mesh  using  a  bilinear  interpolation 
scheme  which  is  identical  to  the  basis  functions  used  in  the  finite  element  model 
(equation  4.9).  This  process  of  refining  the  mesh  and  time  marching  to  steady  state 
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was  continued  three  times  (4  computational  meshes,  total)  until  refinement  of  the  mesh 
did  not  alter  the  model  results.  This  grid  convergence  sequence  resulted  in  the  final 
computational  mesh  shown  in  figure  6.5b.  The  final  mesh  had  1114  nodes  and  1001 
elements.  The  sidewall  transverse  resolution  varied  from  2  elements  at  the  upstream 
end  of  the  modeled  reach  to  3  elements  in  the  vicinity  of  the  transition  beginning  at 
sta  44.31.  The  bottom  transverse  resolution  varied  from  3  elements  at  the  upper  end, 
to  8  elements  at  the  transition,  to  4  elements  at  the  lower  end.  A  transverse-to- 
longitudinal  aspect  ratio  of  1.0  was  used  for  the  elements  within  and  downstream  of 
the  transition.  A  time  step  of  0.1  sec  was  used  for  the  final  grid.  The  primary 
difference  in  the  flow  solutions  provided  by  the  various  grids  was  that  peak  flow 
depths  in  the  oblique  standing  waves  increased  as  the  grid  resolution  increased.  The 
increased  grid  resolution  in  this  region  also  sharpened  the  standing  wave  patterns. 

The  evolution  of  the  flow  depth  at  the  centerline  node  located  at  sta  50.0  is 
presented  in  figure  6.6.  The  length  of  simulation  time  for  each  of  the  four 
computational  meshes  is  also  shown  on  this  figure.  The  time  history  plot  illustrates 
the  time  stepping  to  steady  state  process  at  a  point  in  the  flow  field.  The  depth  at  this 
point  changed  less  than  10'^  ft  after  60  sec. 

Test  Results 

Comparisons  are  made  between  the  water  surface  elevations  measured  in  the 
flume  and  those  obtained  from  the  model  simulations.  Depths  are  reported  relative  to 
the  channel  bottom  as  illustrated  in  figure  6.7.  Flow  depths  in  the  flume  were 
measured  at  small  enough  intervals  so  that  linear  interpolations  of  the  observed  values 
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Figure  6.6.  Time  history  of  depth  for  trapezoidal-to-rectangular  transition  at  centerline 
of  sta  50. 


Figure  6.7.  Definition  sketch  of  depths  plotted  on  depth  contours. 
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provided  an  accurate  map  of  the  water  surface.  The  spatial  intervals  of  observation 
points  varied  from  2.5  ft  longitudinally  and  0.5  ft  transversely  upstream  of  the 
transition  to  as  small  as  0.35  ft  longitudinally  and  0.13  ft  transversely  within  and 
downstream  of  the  transition. 

Contours  of  laboratory  data  and  numerical  model  results  were  plotted  using  the 
FastTABS  postprocessor  (Brigham  Young  University  1993).  The  laboratory  results 
were  placed  on  an  irregular  grid  whose  grid  points  corresponded  to  the  coordinates  at 
which  the  data  were  measured.  First  a  value  at  each  grid  cell  center  is  obtained  using 
a  bilinear  interpolation  of  the  cell’s  nodal  values.  Next,  values  at  the  cell  side 
midpoints  are  obtained  by  linear  interpolation  of  the  nodal  values.  The  cell  is  then 
subdivided  into  a  series  of  triangles  with  vertices  at  the  element  nodes,  at  each 
midpoint  of  the  cell  sides,  and  at  the  cell  center  point.  The  data  are  then  linearly 
interpolated  from  the  values  at  the  vertices.  Finally,  contour  lines  are  drawn  at  a  user 
specified  interval. 

Contours  of  the  observed  and  computed  flow  depths  relative  to  the  channel 
bottom  (h  +  Zb  in  figure  6.7)  are  presented  in  figure  6.8.  Figure  6.9  is  the  computed 
water-surface  mesh.  Water-surface  profiles  along  the  left  wall  and  centerline  are 
provided  in  figure  6.10.  These  results  show  that  the  model  adequately  captures  the 
magnitude  of  the  highest  oblique  standing  waves  but  does  not  accurately  reproduce  the 
wave  phases.  Wave  height  attenuation  through  the  transition  is  abrupt  in  the  shallow 
water  model.  Similar  wave  phase  errors  have  been  observed  by  others  (Jimenez  and 
Chaudhry  1988,  Bhallamudi  and  Chaudhry  1992,  Berger  1992,  Hager  et  al.  1994,  and 
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Berger  and  Stockstill  1994)  in  two-dimensional  shallow  water  applications  to 
supercritical  flow  in  rectangular  channels  and  by  Ippen  (1951)  when  his  analytic 
model  is  compared  with  laboratory  data.  These  wave  phase  errors  are  attributed  to  the 
hydrostatic  pressure  assumption  which  results  in  all  wavelengths  traveling  at  the  speed 


a.  Flume  data 


b.  Model  results 

Figure  6.8.  Depth  contours  (feet  above  channel  bottom)  for  trapezoidal-to-rectangular 


transition. 
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Figure  6.9.  Computed  water-surface  mesh  for  trapezoidal-to-rectangular  transition. 


of  a  long  wave,  c  =  (gh)^^^.  Actually  however,  the  wave  celerity  is  dependent  on  the 
wave  length,  L^,  and  is  given  as  (Lamb  1945) 


bK 

2% 


tanh 


f  ^ 
27rh 


(6.3) 


As  the  depth  to  wavelength  ratio  (h/L^)  becomes  small,  c^  approaches  c. 
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b.  Profile  along  centerline 


Figure  6.10.  Water-surface  profiles  for  trapezoidal-to-rectangular  transition, 
n  =  0.009,  Cu  =  0.1,  =  0.01  ft  (datum  is  10  ft  below  sta  0  invert  elevation) 
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The  standing  waves  patterns  generated  by  supercritical  flow  moving  past  an 
inward  wall  deflection  in  a  rectangular  channel  are  sketched  in  figure  6.11.  The  angle 
between  a  standing  wave  front  and  the  upstream  wall  direction  in  a  rectangular 
channel  is  a  function  of  the  average  approach  velocity  (U)  and  the  wave  celerity 
(Ippen  1951).  For  small  disturbances,  the  actual  wave  angle,  6j^,  is 


(  \ 


5j^  =  sin  * 


(6.4) 


whereas  the  angle  given  by  the  shallow  water  equations  is 


-'sw 


sin 


f 


(6.5) 


STANDING  WAVES 


Figure  6.11.  Plan  view  schematic  of  actual  and  shallow  water  model  standing  wave 


patterns. 
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Computed  and  observed  wavelengths  are  shown  on  the  depth  contour  plots 
presented  in  figure  6.12.  The  measured  wavelength  and  average  flow  depth  in  the 
flume  are  approximately  1.8  ft  and  0.38  ft,  respectively.  This  depth  to  wavelength 
ratio  of  0.21  results  in  an  actual  wave  celerity  of  2.8  ft/sec;  whereas,  the  shallow  water 


a.  Flume  data 


b.  Model  results 

Figure  6.12.  Computed  and  observed  wavelengths  for  trapezoidal-to-rectangular 


transition. 
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wave  celerity  associated  with  a  depth  of  0.38  ft  is  3.5  ft/sec.  At  an  h/L^  of  0.21,  the 
shallow  water  wave  celerity  is  about  25  percent  too  high  due  to  the  inclusion  of  an 
assumed  hydrostatic  pressure  distribution.  Because  the  celerity  computed  by  the 
shallow  water  model  is  too  high,  the  standing  wave  angle  is  too  large  and  the  distance 
between  peak  depths  is  too  short.  Employment  of  higher  order  shallow  water 
equations  which  better  describe  the  fluid  pressure  as  influenced  by  vertical 
accelerations  could  reduce  the  standing  wave  phase  error.  However,  these  Boussinesq 
equations  are  still  limited  to  small  amplitude  waves  (Peregrine  1967)  and  therefore 
will  not  accurately  model  the  wave  patterns  of  three-dimensional  flow  with  large  water 
surface  curvature.  Modeling  considerations  regarding  these  higher  order  equations  are 
further  discussed  in  Chapter  7. 

Velocity  comparisons  are  shown  in  figure  6.13.  Each  reported  measured 
velocity  is  the  average  of  three  point  velocity  measurements  throughout  the  flow  depth 
at  the  station  except  at  stations  near  the  waterline  where  a  single  velocity  at  60  percent 
of  the  depth  was  assumed  to  equal  the  depth  average.  The  three  velocities  were 
measured  at  points  that  were  at  17,  50,  and  83  percent  of  the  depth.  The  velocity 
plots  show  that  the  model  results  compare  well  with  the  observed  velocities. 

The  trapezoidal-to-rectangular  channel  transition  test  illustrates  the  model’s 
ability  to  simulate  oblique  standing  waves.  The  phase  difference  between  the  model 
and  the  real  flow  should  be  of  little  consequence  in  making  engineering  decisions 
regarding  the  maximum  wall  height  required  to  contain  the  flow  in  a  channel  being 
designed.  Because  of  uncertainty  in  the  variations  of  flow  depth  above  temporal  mean 
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LONGITUDINAL  VELOCITY,  ft/sec 


LONGITUDINAL  VELOCITY,  ft/sec 


a.  sta  30.0 


b.  sta  40.0 


Figure  6.13.  Computed  and  observed  depth-averaged  longitudinal  velocities  for 
trapezoidal-to-rectangular  transition  (left  and  right  referenced  to  looking  downstream). 


averages,  constructed  wall  heights  are  a  maximum  in  the  vicinity  of  transitions  and  are 
only  gradually  decreased  toward  stations  where  the  flow  is  near  uniform,  thus 
providing  adequate  freeboard  throughout  the  region  of  large  standing  waves. 

Parameter  Sensitivity  Tests 

Additional  simulations  of  the  trapezoidal-to-rectangular  transition  were 
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conducted  to  evaluate  the  model’s  sensitivity  to  the  values  of  specified  parameters.  In 
each  of  these  tests,  one  parameter  was  changed  from  the  previously  reported  values 
(table  6.2)  and  the  model  was  time  stepped  to  the  steady  state  solution  using  the 
results  shown  in  figure  6.8  as  the  initial  conditions.  The  parameters  evaluated  were 
the  Manning’s  roughness  coefficient  (n),  the  turbulent  kinematic  viscosity  coefficient 
(Cj,),  and  the  constant  depth  at  the  moving  boundary  nodes  (fi). 

Test  1;  Manning’s  n  Value 

The  Manning’s  n  value  was  changed  from  0.009  to  0.010.  The  steady  state 
solution  for  n  =  0.010  is  shown  as  water-surface  profiles  in  figure  6.14.  The  larger  n 
value  resulted  in  a  hydraulic  jump  forming  near  sta  48,  just  upstream  of  the  lower  end 
of  the  transition.  The  energy  losses  along  the  transition  approach  with  the  increased 
bed  roughness  resulted  in  insufficient  specific  energy  to  pass  the  flow  through  the 
transition  as  supercritical.  Examination  of  the  water  surface  elevations  upstream  of  the 
transition  indicates  that  the  computed  flow  depth  is  larger  than  that  measured  in  the 
flume;  whereas,  the  computed  flow  depth  for  a  n  =  0.009  is  very  close  to  the  observed 
depth  (figure  6.10).  These  types  of  laboratory  flumes  have  been  documented  at  the 
Waterways  Experiment  Station  as  having  a  roughness  equivalent  to  a  n  =  0.009; 
however,  data  reduction  of  these  tests  were  based  on  unidirectional  flow  using  the 
hydraulic  radius  as  the  flow  area  divided  by  the  wetted  perimeter  rather  than 
employing  the  wide  channel  assumption  that  the  flow  depth  approximates  the 
hydraulic  radius  as  is  used  in  the  two-dimensional  model.  The  Test  1  results  indicate 
that  n  =  0.009  is  appropriate  for  two-dimensional  modeling  of  the  epoxy  painted 
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a.  Profile  along  left  wall 


b.  Profile  along  centerline 

Figure  6.14.  Water-surface  profiles  for  trapezoidal-to-rectangular  transition, 
n  =  0.010,  Cb  =  0.1,  ft  =  0.01  ft  (datum  is  10  ft  below  sta  0  invert  elevation). 
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surface  of  the  channel  shapes  and  sizes  constructed  for  the  flume  tests  used  here. 

Test  2:  Turbulent  Kinematic  Viscosity  Coefficient 

The  turbulent  kinematic  viscosity  coefficient  (C|,)  used  in  4.5  was  increased  by 
an  order  of  magnitude  to  test  the  model  sensitivity  to  selection  of  this  parameter.  A 
Cj,  value  of  1.0  was  chosen  resulting  in  a  turbulent  kinematic  viscosity  (V()  of  0.23 
ft^/sec.  Initial  conditions  for  this  test  were  the  steady  state  solution  shown  in  figure 
6.8.  The  order  of  magnitude  increase  in  did  not  have  a  significant  effect  on  the 
steady  state  solution  for  this  particular  channel  and  flow  configuration.  Water-surface 
profiles  resulting  from  a  Vj  of  0.23  ft^/sec  are  shown  in  figure  6.15.  The  increased 
viscosity  did  not  affect  the  amplitude  of  the  first  standing  wave  at  the  channel 
centerline  (figure  6.15b)  which  compares  well  with  the  laboratory  data.  However,  use 
of  the  increased  viscosity  dampened  the  computed  reflected  waves  downstream  of  the 
first  waves.  The  flow  depth  approaching  the  transition  seems  to  be  modeled  correctly 
using  this  large  viscosity,  yet  the  relatively  smooth  water  surface  downstream  of  the 
first  wave  is  unrealistic.  These  test  results  show  that  a  Cj,  value  of  0.1  provides  much 
better  model  results  than  a  C|,  value  of  1.0.  The  viscosity  associated  with  the  large 
value  of  Cj,  over  damped  the  standing  waves  following  the  initial  wave. 

Test  3:  Flow  Depth  at  the  Moving  Boundary 

The  final  sensitivity  test  evaluated  the  model  response  to  a  smaller  flow  depth 
at  the  moving  boundary  nodes.  A  depth  of  fi  =  0.005  ft  was  selected  for  this  test. 
This  constant  depth  is  one  half  the  value  used  for  the  previously  discussed  tests  (figure 
6.8).  Water-surface  profiles  of  the  steady  state  solution  with  equal  to  0.005  ft  are 
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ile  along  left  wall 


b.  Profile  along  centerline 

Figure  6.15.  Water-surface  profiles  for  trapezoidal-to-rectangular  transition, 
n  =  0.009,  C|,  =  1.0,  fi  =  0.01  ft  (datum  is  10  ft  below  sta  0  invert  elevation). 


104 


shown  in  figure  6.16.  As  expected,  the  water  surface  elevations  are  practically 
identical  to  those  computed  using  an  fi  value  of  0.01  ft. 

At  this  point  in  the  testing  process,  the  model’s  sensitivity  to  selection  of  initial 
conditions  was  investigated.  The  model’s  sensitivity  to  each  of  the  three  parameters 
was  reevaluated  using  a  grid  similar  in  resolution  to  that  shown  in  figure  6.5b,  but 
rather  than  using  the  flow  conditions  shown  in  figure  6.8  as  the  initial  conditions, 
uniform  flow  conditions  (normal  depth  and  velocity)  were  used  as  the  initial 
conditions.  The  steady  state  solution  resulting  from  these  initial  conditions  were  for 
all  practical  purposes  indistinguishable  from  the  results  reported  above  in  the  Test  1 
through  Test  3  sections. 

These  results  indicate  that  for  the  discharge  and  geometric  design  tested,  the 
model  is  insensitive  to  the  choice  of  flow  depth  at  the  moving  boundary  within  the 
range  tested  and  that  it  is  only  mildly  sensitive  to  the  choice  of  turbulent  kinematic 
viscosity  coefficient.  However,  for  the  configuration  tested,  the  model  was  found  to 
be  extremely  sensitive  to  the  selection  of  Manning’s  n.  Therefore,  in  practice  the 
hydraulic  design  process  should  include  the  evaluation  of  flow  conditions  resulting 
from  a  range  of  reasonable  roughness  coefficients  such  that  an  envelope  of  possible 
hydraulic  conditions  can  be  established.  In  a  hydraulic  design  situation,  without  the 
benefit  of  laboratory  observations,  the  model  indicated  that  the  1  lateral  on  6 
longitudinal  transition  is  too  abrupt  and  that  it  may  result  in  a  choked  flow  condition. 

When  the  model  results  are  compared  to  the  observed  flume  data,  it  is  apparent 
that  initially  chosen  values  for  these  parameters  (table  6.2)  were  appropriate  for 
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b.  Profile  along  centerline 

Figure  6.16.  Water-surface  profiles  for  trapezoidal-to-rectangular  transition, 
n  =  0.009,  Q  =  0.1,  fi  =  0.005  (datum  is  10  ft  below  sta  0  invert  elevation). 
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evaluation  of  the  hydraulic  performance  of  the  channel.  Therefore,  these  values  were 
used  in  aU  subsequent  model  testing. 

Trapezoidal  Channel  with  a  Horizontal  Curve 

The  flume,  part  of  a  site  specific  1: 16-scale  physical  model  study  of  a  flood 
control  channel  initiated  for  the  U.S.  Army  Engineer,  Los  Angeles  District,  was  a 
trapezoidal  channel  having  a  horizontal  curve.  A  portion  of  the  modeled  reach 
consisted  of  a  hydraulically  steep,  trapezoidal  flume  having  a  long  horizontal  curve. 
The  bend  provided  a  good  test  case  since  the  direction  of  nodal  displacement  in  the 
numerical  simulation  would  vary  along  the  length  of  the  bend  and  the  flow  boundary 
on  the  inside  of  the  bend  would  have  a  different  elevation  from  that  on  the  outside  due 
to  the  superelevation  of  the  water  surface. 

The  portion  of  the  facility  tested  for  this  study  was  nominally  52.84  ft  long  on 
centerline,  with  wall  heights  ranging  from  0.31  ft  to  0.44  ft  above  the  invert,  with  base 
widths  varying  from  0.938  ft  to  0.625  ft  and  side  slopes  varying  from  1  vertical  on  3 
horizontal  to  1  vertical  on  2  horizontal.  The  channel  slope  was  0.0073.  The  side 
slopes  in  the  8.0-ft-long  straight  reach  (sta  64.741  to  sta  56.741)  and  in  the  5.313-ft- 
long  transition  section  (sta  56.741  to  sta  51.428)  were  covered  with  crushed  rock 
having  an  average  diameter  of  3/8  inches.  The  flow  width  within  the  straight  reach 
was  bounded  to  2.813  ft  by  vertical  walls  at  the  top  of  the  rock  slopes.  All  other 
boundaries  were  constructed  of  plastic  coated  plywood  covered  with  rolled  epoxy 
paint.  The  channel  bend  was  composed  of  a  compound  curve  in  which  spiral  curves 
were  used  to  transition  between  the  upstream  and  downstream  tangents  to  a  circular 
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curve.  The  circular  curve  was  of  radius  15.625  ft.  A  plan  view  of  the  flume  is  shown 
in  figure  6.17  and  photographs  of  the  flume  are  shown  in  figure  6.18. 


Figure  6.17.  Plan  view  and  sections  of  trapezoidal  horizontal  curve. 


a.  Looking  downstream 

Figure  6.18.  Dry  bed  view  of  trapezoidal  horizontal  curve. 
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Water  was  supplied  from  a  pump  to  an  8-inch  diameter  inflow  pipe.  The 
flume  discharge  was  returned  to  the  sump  through  a  drain  line.  The  rate  of  flow 
entering  the  flume  was  measured  using  a  6-inch  by  3.6-inch  venturi  meter  installed  in 
the  inflow  pipe.  Point  gages  were  used  to  measure  flow  depths  and  a  propeller  type 
probe  was  used  to  measure  velocities.  Flow  was  introduced  in  the  flume  by  means  of 
a  headbox.  The  flume  inflow  was  subcritical  and  therefore  a  sluice  gate  was  not 
required.  The  subcritical  inflow  was  established  by  downstream  control  within  the 
rock  lined  transition  where  the  flow  accelerated  through  critical.  Upstream  of  the 
transition  the  flow  was  subcritical  and  supercritical  flow  existed  downstream  of  the 
transition.  The  outflow  downstream  of  the  transition  was  supercritical  and  therefore 
no  tailwater  level  was  set. 

Test  Conditions 

Comparisons  were  made  of  model  results  with  lab  data  obtained  in  the  flume 
with  a  discharge  of  1.22  ft^/sec.  The  flow  conditions  in  the  laboratory  flume  are 
shown  in  figure  6.19.  The  numerical  model  used  a  Manning’s  n  of  0.009  for  the 
epoxy  painted  plywood  area.  The  n  value  associated  with  the  crushed  rock  placed  on 
the  side  slopes  at  the  upper  end  was  estimated  to  be  0.021  using  the  following  form  of 
Strickler’s  equation  (US  Army,  Office,  Chief  of  Engineers  1991): 

n  =  0.036  (6.6) 

where  is  the  size  of  which  90  percent  of  the  sample  is  finer  and  is  given  in  feet. 
The  value  for  D^q  was  estimated  assuming  a  linear  variation  in  gradation  from  the 
Dioo  to  the  D50  values.  The  roughness  parameters  were  chosen  prior  to  the  numerical 


Figure  6.19.  Flow  conditions  in  trapezoidal  horizontal  curve  (looking  upstream). 

model  simulations.  The  values  were  never  changed  (or  adjusted)  to  make  the  model 
results  agree  with  the  measured  data;  the  model  was  not  calibrated  to  reproduce  the 
flume  data. 

As  with  the  first  test,  the  initial  top  width  of  the  flow  was  based  on  calculated 
normal  depths.  Normal  depths  were  computed  assuming  uniform  straight  channels 
having  cross  sectional  shapes  and  slopes  identical  to  the  curved  flume.  Two  depths 
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were  computed,  one  for  the  upper  end  using  an  n  value  of  0.021  for  a  channel  having 
a  base  width  of  0.938  ft  and  1  vertical  on  3  horizontal  side  slopes  and  one  using  an  n 
value  of  0.009  for  a  channel  having  a  base  width  of  0.625  ft  and  1  vertical  on  2 
horizontal  side  slopes.  Both  uniform  depths  were  computed  using  the  channel  slope  of 
0.0073.  Given  a  flow  depth,  the  channel  base  width,  and  side  slopes,  the  top  width 
was  computed  for  the  initial  grid.  The  normal  depth  for  the  upper  end  was  0.31  ft  and 
normal  velocity  was  2.08  ft/sec  (Froude  number  of  0.81).  The  normal  depth  for  the 
section  having  1  vertical  on  2  horizontal  side  slopes  was  0.25  ft  and  the  normal 
velocity  was  4.25  ft/sec  (Froude  number  of  1.79).  The  initial  mesh,  shown  in  figure 
6.20a,  had  629  nodes  and  541  elements  and  was  considered  coarse,  but  adequate  for 
startup.  Lateral  resolution  of  the  channel  bottom  was  3  elements  for  the  upstream 
portion  (bottom  width  =  0.938  ft)  and  2  elements  for  the  remainder  of  the  channel 
(bottom  width  =  0.625  ft).  Two  elements  were  used  on  each  side  slope  of  the  initial 
mesh. 

In  the  same  manner  as  the  previous  trapezoidal-to-rectangular  transition  was 
modeled,  the  upstream  and  downstream  ends  of  the  mesh  were  composed  of  additional 
points  representing  a  headbox"  in  addition  to  points  added  at  the  downstream  end 
representing  a  "tailbox".  These  additional  points  were  used  to  simplify  the  boundary 
condition  input.  A  rectangular  section  at  the  inflow  boundary  simplifies  the 
specification  of  inflow  discharge.  The  unit  discharge  is  merely  the  total  discharge 
divided  by  the  width  of  the  rectangular  inflow  boundary;  whereas,  the  unit  discharge  at 
a  trapezoidal  inflow  section  varies  with  depth  assuming  uniform  velocity  across  the 
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Figure  6.20.  Numerical  model  computational  meshes  for  trapezoidal  horizontal  curve. 

inflow  boundary.  The  use  of  a  rectangular  section  at  the  outflow  boundary  results  in 
fixed  nodes  at  the  waterline  at  the  outflow  boundary.  Because  a  moving  waterline 
boundary  at  the  outflow  boundary  produces  an  additional  source  of  nonlinearity  in  the 
solution,  a  trapezoidal  section  at  the  outflow  boundary  is  not  desirable.  This  is  the 
reason  additional  elements  representing  a  tailbox  were  used  at  the  outflow  boundary 
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which  simplified  the  model’s  adherence  to  specified  outflow  boundary  conditions. 
These  "boxes"  in  which  the  boundary  nodes  were  fixed  produced  a  gentle  transition 
from  rectangular  to  trapezoidal  shape  at  the  upstream  end  and  trapezoidal  to 
rectangular  shape  at  the  downstream  end. 

The  inflow  boundary  flow  was  subcritical  and  required  the  specification  of  p 
and  q  at  each  inflow  node.  Values  of  p  and  q  were  determined  as  the  total  discharge 
(1.22  ftVsec),  applied  normal  to  the  inflow  boundary,  divided  by  the  width  of  the 
inflow  boundary  (2.758  ft).  The  outflow  boundary  was  supercritical  and  therefore,  no 
boundary  conditions  were  specified.  The  upstream  and  downstream  boundary 
conditions  are  listed  in  table  6.3  where  the  values  of  p  and  q  are  valid  for  the  x  and  y 
coordinate  system  orientation  shown  in  figure  6.20.  The  coordinate  system  orientation 
coincided  with  the  North  and  East  directions  of  the  prototype  project  plans  from  which 
the  flume  was  constructed.  The  model  parameters  used  for  this  simulation  are 
presented  in  table  6.4.  The  specified  boundary  depth,  1i  =  0.01  ft,  was  less  than  5 
percent  of  the  depth  at  the  channel  center  in  each  of  the  two  reaches  thus  producing  a 
cross  sectional  area  error  of  roughly  0.06  percent.  This  specified  boundary  depth  of 
0.01  ft  applied  to  normal  flow  conditions  in  section  A-A  (figure  6.17)  resulted  in  an 
initial  grid  width  of  2.258  ft  for  this  section.  An  initial  grid  width  of  1.705  ft  for 

section  B-B  was  computed  using  normal  depth  (0.28  ft)  and  a  boundary  depth  of  0.01 
ft. 

Initially  0.1  sec  time  steps  were  used  for  the  simulation  for  a  time  of  15  sec. 
This  time  step  resulted  in  a  Courant  number  of  1.0.  The  time  steps  were  then 
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Table  6.3:  Curved  channel,  trapezoidal  cross  section,  boundary  conditions. 


Boundary 

Flow  Regime 

Boundary  Condition 

Value 

Inflow 

Subcritical 

h 

Not  Specified 

P 

-0.1999  tf/sec 

q 

-0.3947  ftVsec 

Outflow 

Supercritical 

h 

Not  Specified 

P 

Not  Specified 

q 

Not  Specified 

Table  6.4:  Curved  channel,  trapezoidal  cross  section,  model  parameters  (Uj  is 
Manning’s  coefficient  for  the  smooth  plywood  and  nj  is  for  the  crushed  rock). 


Condition 

Value 

a 

1.0 

p 

0.1 

Hi,  nj 

0.009,  0.021 

Cb,  V, 

0.1,  0.015  ftVsec 

0.01  ft 
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increased  to  0.25  sec  and  simulations  were  conducted  for  an  additional  235  seconds 
resulting  in  a  total  time  of  250  sec. 

The  computed  flow  variables  and  the  boundary  location  at  time  250  sec  were 
equal  to  those  resulting  from  245  sec  of  simulation  and  therefore,  this  solution  was 
taken  to  represent  the  steady  state.  After  the  steady  state  solution  was  obtained,  the 
initial  mesh  was  refined.  Initial  conditions  for  the  refined  mesh  were  obtained  from 
the  coarser  mesh  results  using  a  bilinear  interpolation  scheme.  Simulations  were  then 
run  until  steady  state  was  obtained  with  the  refined  mesh.  This  process  was  continued 
until  mesh  convergence  was  reached,  this  being  the  point  at  which  the  solution  did  not 
change  with  additional  mesh  refinement.  In  this  process  of  mesh  convergence,  a  time 
step  of  0.25  sec  was  used  for  a  total  simulation  time  of  385  sec.  The  resulting  mesh 
is  shown  in  figure  6.20b.  This  final  mesh  had  1450  nodes  and  1303  elements.  The 
primary  difference  in  the  coarse  and  resolved  grids  was  that  the  finer  grid  provided 
slightly  larger  peak  superelevations  in  the  water  surface  resulting  from  runup  and 
drawdown  on  the  outside  and  inside  of  the  bend,  respectively. 

Test  Results 

A  plot  of  the  computed  and  observed  side  boundaries  is  provided  in  figure 
6.21.  The  model  accurately  solved  the  waterlines  through  the  transition  where  the 
flow  accelerated  from  subcritical  to  supercritical.  Contours  of  the  flow  depth  relative 
to  the  channel  bottom  are  shown  for  the  laboratory  data  and  simulation  results  in 
figure  6.22.  Longitudinal  spacing  of  flow  depth  measurement  varied  from  1  ft  to 
0.27  ft  while  transverse  spacing  varied  from  0.31  ft  near  the  channel  centerline  to  as 
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Figure  6.21.  Computed  and  observed  side  boundaries  for  trapezoidal  horizontal  curve. 


small  as  0.02  ft  near  the  waterline.  Depth  contours  of  the  laboratory  data  (figure 
6.22a)  show  that  the  entrance  conditions  caused  flow  disturbances  at  the  upper  end  of 
the  flume.  The  poor  entrance  conditions  were  of  little  importance  since  the  primary 
concern  of  this  study  is  the  flow  conditions  within  the  channel  curve.  It  is  interesting 
to  note  the  computed  water-surface  oscillations  downstream  of  the  bend  as  the  flow 
adjusts  to  uniformity.  Unfortunately  the  tangent  reach  in  the  flume  was  not  long 
enough  to  produce  these  disturbances  caused  by  the  channel  bend.  Cross  sectional 
plots  of  measured  and  computed  water  surface  elevations  near  the  midpoint  of  the 
curve,  at  the  downstream  end  of  the  simple  curve,  and  at  the  downstream  end  of  the 
spiral  curve  are  provided  in  figure  6.23.  These  results  show  that  the  superelevation  in 
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Figure  6.23.  Cross  section  view  of  flow  depth  for  trapezoidal  horizontal  curve 
(looking  downstream). 
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the  bend  is  adequately  captured.  Figure  6.22  shows  that  the  model  adequately 
reproduces  the  amplitude  of  wave  oscillations. 

Depth-averaged  longitudinal  velocities  at  three  stations  are  plotted  in  figure 
6.24.  The  depth-averaged  velocities  were  generally  orientated  in  directions  parallel  to 
the  channel  centerline.  Depth-averaged  velocities  were  obtained  in  the  laboratory  by 
measuring  three  point  velocities  (at  17,  50,  and  83  percent  of  the  depth)  throughout  the 
flow  column  except  at  stations  near  the  waterline  where  the  velocity  measured  at  60 
percent  of  depth  was  assumed  to  represent  the  depth  average.  The  entrance  conditions 
probably  produced  the  unsymmetrical  velocity  distribution  at  sta  51.43.  The  depth 
contours  and  velocity  plots  indicate  that  the  computed  depths  at  sta  51.43  were  too 
large  and  the  velocities  too  low.  However,  the  differences  are  small  and  suggest  that 
the  Manning’s  coefficient  used  for  the  crushed  rock  was  perhaps  too  large.  The  lateral 
velocity  distribution  at  sta  43.62  (beginning  of  the  circular  curve)  is  uniform  in  the 
channel  center,  but  the  largest  velocities  have  migrated  toward  the  outside  of  the  bend 
at  sta  30.90  (approximately  the  midpoint  of  the  curve). 

The  depth-averaged  model  is  not  capable  of  simulating  the  migration  of 
maximum  streamwise  velocity  toward  the  outside  of  the  bend  because  the  so-called 
secondary  flow  that  occurs  in  channel  bends  is  a  three  dimensional  phenomenon  and  is 
not  addressed  using  the  shallow  water  equations.  Because  of  this  limitation, 
engineering  decisions  such  as  riprap  size  selection  for  bank  protection  in  trapezoidal 
channel  bends  should  not  be  based  solely  on  the  results  of  depth-averaged  models. 
These  models  can  not  provide  the  location  and  magnitude  of  the  maximum  channel 
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Figure  6.24.  Computed  and  observed  depth-averaged  velocities  for  trapezoidal 
horizontal  curve  (left  and  right  directions  are  referenced  to  looking  downstream). 
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velocity  in  trapezoidal  channel  bends.  However,  the  model  does  capture  the 
superelevation  in  the  water  surface  within  the  channel  bend  which  is  most  important  in 
flood  control  channel  design.  Although  the  model  did  not  accurately  reproduce  the 
velocities  near  the  midpoint  of  the  curve  (figure  6.24c)  the  water  surface  elevations  in 
this  vicinity  were  accurately  computed  (figure  6.23a).  Water  surface  elevations  are  of 
primary  concern  in  designing  high-velocity  channels  and  the  model  has  been  shown  to 
predict  these  satisfactorily. 

Rectangular-to-Trapezoidal  Transition 
Additional  model  testing  was  conducted  to  evaluate  the  model’s  ability  to 
capture  oblique  standing  waves  in  trapezoidal  channels.  Data  were  obtained  in  the 
tilting  open  channel  flume  previously  described  (figure  6.1)  having  a  base  width  of  2.0 
ft  and  side  slopes  of  1  vertical  on  2.25  horizontal.  A  rectangular-to-trapezoidal 
channel  transition  was  constructed  by  placing  vertical  walls  at  the  toe  of  the  side 
slopes  along  the  first  30.73  ft  of  the  trapezoidal  flume.  These  vertical  walls  were  then 
continued  up  the  side  slopes  for  a  distance  of  13.5  ft  resulting  in  a  1  lateral  to  6 
longitudinal  transition  on  each  sidewall  (figure  6.25).  The  resulting  flume  cross 
sections  were  of  rectangular  shape  from  the  upstream  end  at  sta  0.0  to  sta  30.73  and  of 
trapezoidal  shape  from  sta  44.23  to  the  flume  end  at  sta  80.0.  A  dry  bed  photograph 
of  the  rectangular-to-trapezoidal  transition  is  provided  in  figure  6.26. 

Test  Conditions 

The  flow  conditions  documented  in  the  flume  experiments  and  subsequently 
simulated  using  the  numerical  model  consisted  of  a  discharge  of  3.6  ftVsec  on  a  slope 
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Figure  6.25.  Plan  view  and  sections  of  rectangular-to-trapezoidal  transition. 

of  0.012.  Given  uniform  conditions,  this  channel  configuration  results  in  a  computed 
Froude  number  of  2.2  in  the  transition  approach  where  the  channel  was  a  2-ft  wide 
rectangular  section.  Flow  conditions  within  and  downstream  of  the  transition  are 
shown  in  figure  6.27.  As  the  flow  flared  out  at  the  upstream  end  of  the  transition,  the 
flow  on  each  side  met  the  steep  side  slopes.  This  produced  a  rise  in  the  water  surface 
within  the  portion  of  the  transition  located  over  the  sloping  sidewalls.  The  flow  depth 
near  the  channel  centerline  decreased  beginning  at  the  upstream  end  of  the  transition 
due  to  the  expansion  of  the  channel.  The  water  runup  on  the  side  slopes  and  the 
depression  along  the  centerline  were  followed  by  a  depression  along  the  side  slopes 
and  a  rapid  increase  along  the  centerline  as  the  standing  waves  generated  at  the  side 
slopes  intersected.  These  oblique  standing  waves  resulted  in  varying  waterlmes  along 

the  left  and  right  side  slopes. 


Figure  6.26.  Dry  bed  view  of  rectangular-to-trapezoidal  transition  (looking 
downstream). 
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In  the  same  manner  as  the  first  two  tests,  the  initial  flow  conditions  and  initial 
grid  width  were  based  on  normal  flow  conditions  assuming  a  Manning’s  n  value  of 
0.009.  Normal  depth  (0.28  ft)  and  velocity  (6.5  ft/sec)  were  computed  for  the 
rectangular  section  and  used  as  initial  conditions  for  the  upstream  reach  (sta  0.0  to  sta 
30.73).  Also,  normal  depth  (0.24  ft)  and  velocity  (6.0  ft/sec)  for  the  trapezoidal 
section  were  used  as  initial  conditions  for  the  downstream  reach  (sta  34.835  to  sta 
80.0).  The  constant  depth  at  the  moving  boundaries  within  the  trapezoidal  section,  fi, 
was  set  to  0.01  ft  which  was  about  4  percent  of  normal  depth.  The  normal  depth  and 
were  used  to  compute  the  initial  mesh  width  (3.035  ft)  within  the  trapezoidal  reach. 
Average  values  of  the  flow  conditions  in  the  upper  and  lower  reaches  were  applied  as 
initial  conditions  within  the  transition.  The  initial  computational  mesh  (figure  6.28a) 
had  846  nodes  and  759  elements.  The  lateral  resolution  included  3  elements  at  the 
upstream  end  of  the  rectangular  reach.  The  mesh  contained  5  elements  across  the 
rectangular  channel  immediately  upstream  of  the  transition  beginning  at  sta  26.34, 
across  the  base  of  the  transition,  and  across  the  base  of  the  trapezoidal  channel 
downstream  of  the  transition.  Three  elements  were  used  across  each  of  the  side  slopes 
of  the  trapezoidal  reach  from  the  end  of  the  transition  to  sta  72.6.  At  the  lower  end  of 
the  flume,  the  mesh  had  2  elements  across  each  side  slope  and  3  elements  across  the 
base.  Additional  nodes  and  elements  were  added  to  the  downstream  end  representing 
a  4.41 -ft  long  gentle  transition  from  a  trapezoidal  to  a  rectangular  section.  The 
outflow  boundary  elements  represented  a  rectangular  "tailbox";  the  side  wall  nodes  at 
the  outflow  boundary  were  fixed. 
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Figure  6.28.  Numerical  model  computational  meshes  for  rectangular-to-trapezoidal 
transition. 


Flow  at  the  inflow  boundary  was  supercritical  requiring  input  of  p,  q,  and  h. 
Since  the  computational  mesh  was  aligned  with  the  x  axis  (figure  6.28),  values  of  q 
for  each  inflow  node  were  set  to  zero.  The  p  values  at  these  nodes  were  set  equal  to 
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the  flow  rate  per  unit  width  of  channel  (total  discharge,  3.6  ftVsec,  divided  by  channel 
width,  2  ft).  As  mentioned  above,  depth  at  the  inflow  nodes  was  set  to  normal  depth. 
The  outflow  was  supercritical  and  therefore  no  boundary  conditions  were  specified. 
The  boundary  conditions  employed  in  the  rectangular-to-trapezoidal  transition  are 
shown  in  table  6.5  and  the  model  parameters  used  are  listed  in  table  6.6. 

Boundary  nodes  at  the  waterline/transition  wall  intersection  on  each  side  of  the 
channel  were  allowed  to  move  in  the  direction  parallel  to  the  vertical  transition  walls. 
As  the  water  surface  elevation  at  this  point  increased  or  decreased,  the  intersection 
point  moved  downstream  or  upstream,  respectively.  All  other  side  slope  boundary 
nodes  were  specified  to  move  in  the  direction  of  maximum  side  slope.  Each  of  the 
nodes  within  the  rectangular  reach  was  fixed. 

Table  6.5:  Rectangular-to-trapezoidal  transition,  boundary  conditions. 


Boundary 

Flow  Regime 

Boundary  Condition 

Value 

Inflow 

Supercritical 

h 

0.280  ft 

P 

1.80  ft^/sec 

q 

0.0 

Outflow 

Supercritical 

h 

Not  Specified 

P 

Not  Specified 

q 

Not  Specified 
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Table  6.6:  Rectangular-to-trapezoidal  transition,  model  parameters. 


Condition 

Value 

a 

1.0 

p 

0.1 

n 

0.009 

Cb-  Vt 

0.1,  0.023  ftVsec 

ft 

0.01  ft 

With  these  initial  conditions,  boundary  conditions,  and  nodal  movement 
constraints,  4.0  sec  of  simulation  were  run  using  a  time  step  of  0.05  sec  which  resulted 
in  a  Courant  number  of  1.3.  A  time  step  of  0.1  sec  was  then  used  for  a  total 
simulation  time  of  19.0  sec.  At  this  time  the  mesh  was  refined  such  that  the  grid 
density  upstream,  through,  and  downstream  of  the  transition  was  increased. 

The  final  grid  had  1406  nodes  and  1297  elements.  The  mesh  included  7 
elements  across  the  rectangular  reach  immediately  upstream  of  the  transition  and 
across  the  base  within  the  transition,  resulting  in  an  element  longitudinal-to-transverse 
aspect  ratio  of  approximately  1.5.  Seven  elements  were  used  across  the  base  and  4 
elements  across  each  side  slope  from  the  transition  to  sta  51.3.  The  grid  was  not 
further  refined  below  sta  51.3.  Time  steps  varying  from  0.01  sec  to  0.2  sec  were  used 
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to  simulate  a  total  time  of  28.6  sec,  the  results  of  which  were  deemed  steady  state. 
The  final  grid  corresponding  to  the  steady  state  solution,  is  shown  in  figure  6.28b. 

The  aspect  ratios  of  the  elements  located  on  the  side  slope  varied  along  the  length  of 
the  channel  in  conjunction  with  the  waterline  variance  (transition  detail  in  figure 
6.28b).  In  the  vicinity  of  maximum  flow  width,  the  side  slope  elements  had  aspect 
ratios  (longitudinal  to  lateral)  of  approximately  3;  whereas  in  regions  of  minimum 
flow  width  the  element  aspect  ratios  were  about  5.  These  aspect  ratios  are  appropriate 

for  modeling  high-velocity  channel  fiow  in  rectangular  channels  (Stockstill  and  Berger 
1994). 

Test  Results 

Computed  water  surface  elevations  are  compared  with  those  measured  in  the 
laboratory  flume  by  examination  of  the  computed  and  observed  depth  contours  shown 
in  figure  6.29.  These  depths  are  relative  to  the  channel  bottom  (see  figure  6.7).  The 
observed  values  are  interpolations  from  discrete  measurements  obtained  at  a  1.0  ft 
longitudinal  spacing  upstream  of  the  transition  to  as  small  as  0.5  ft  in  the  vicinity  of 
sta  39.  Lateral  distances  between  depth  measurements  varied  from  1.0  ft  upstream  of 
the  transition  (at  each  wall  and  on  centerline)  to  0.07  ft  in  the  vicinity  of  the  standing 
waves.  The  model  results  are  further  illustrated  in  the  water-surface  mesh  presented  in 
figure  6.30. 

Evaluation  of  the  observed  and  computed  depth  contours  (figure  6.29), 
waterlines  (figure  6.31),  and  water-surface  profiles  (figure  6.32)  leads  to  the 
conclusion  that  in  this  flow  situation  also  there  is  a  phase  difference  in  the  computed 


131 


a.  Flume  data 


b.  Model  results 

Figure  6.29.  Depth  contours  (feet  above  channel  bottom)  for  rectangular-to-trapezoidal 
transition. 


and  observed  wave  patterns.  The  computed  location  of  the  first  wave  peak  is  slightly 
downstream  of  that  measured  in  the  flume.  This  is  opposite  from  the  test  results  for 
the  trapezoidal-to-rectangular  transition  where  the  standing  waves  were  contained  in 
the  rectangular  portion  of  the  channel.  This  error  is  probably  due  to  the  violation  of  a 
combination  of  model  assumptions  that  the  slope  is  geometrically  mild,  that  the 
pressure  is  hydrostatic,  and  that  the  velocity  is  uniformly  distributed  throughout  the 


fluid  column. 
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Figure  6.30.  Computed  water-surface  mesh  for  rectangular-to-trapezoidal  transition. 
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Figure  6.31.  Computed  and  observed  side  boundaries  for  rectangular-to-trapezoidal 
transition. 


The  computed  standing  wave  phase  error  may  be  attributed  in  part  to  the  depth 
averaging  of  velocity  assuming  a  uniform  vertical  distribution  of  horizontal  velocity. 
The  currents  in  the  rectangular  reach  of  the  flume  are  larger  near  the  surface  than  near 
the  bed  due  to  boundary  drag.  A  short  distance  downstream  of  the  vertical  wall  break, 
the  flow  near  the  surface  is  directed  approximately  parallel  with  the  vertical  wall  while 
the  currents  near  the  slope  toe  are  still  generally  in  the  longitudinal  direction  resulting 
in  secondary  currents  in  this  vicinity.  There  is  a  significant  vertical  variation  in 
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a.  Profile  along  left  wall 


b.  Profile  along  centerline 

Figure  6.32.  Water-surface  profiles  for  rectangular-to-trapezoidal  transition  (datum  is 
10  ft  below  sta  0  invert  elevation). 
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horizontal  velocity  throughout  the  flow  near  the  transition  wall,  but  the  depth-averaged 
model  does  not  account  for  this  vertical  variation  in  velocity.  This  discrepancy  could 
contribute  to  the  standing  wave  phase  differences  between  the  depth-averaged  model 
and  the  flume.  However,  it  is  difficult  to  quantify  how  these  velocity  errors  impact 
the  computed  water  surface  elevation.  One  can  speculate  that  the  depth-averaged 
velocity  given  by  the  model,  assuming  no  vertical  variation  in  horizontal  velocity,  is 
faster  than  that  in  the  real  system  (flume)  near  the  transition  wall.  This  higher 
velocity  would  coincide  with  a  shallower  depth,  as  shown  in  the  water-surface  profile 
along  the  wall  (figure  6.32a).  Perhaps  the  computed  momentum  of  the  flow  moving 
up  the  side  slope  is  excessive  resulting  in  the  peak  runup  occurring  too  far 
downstream. 

Another  possible  explanation  of  model  and  flume  differences  is  related  to  the 
mild-slope  assumption  employed  in  the  model  equations  because  as  the  flow  expands 
and  contracts,  the  bed  slope  for  flow  near  the  waterline  is  significantly  steeper  than  the 
channel’s  longitudinal  slope.  A  qualitative  assessment  of  the  error  introduced  by  the 
mild-slope  assumption  can  be  made  by  considering  a  steep  channel  in  which  the  flow 
is  nearly  parallel  to  the  bed  (see  figure  6.33).  Under  these  circumstances,  the  steep- 
flow  equations  are  a  good  description  of  the  flow  since  accelerations  normal  to  the  bed 
are  nearly  zero,  although  vertical  accelerations  may  be  significant.  Insight  into 
solutions  given  by  the  mild-slope  form  of  the  shallow  water  equations  when  applied  to 
steep  channels  can  be  obtained  by  substitution  of  the  mild-slope  variables  into  the 
steep-slope  equations  while  assuming  the  steep-slope  equations  perfectly  describe  the 


The  steady  one-dimensional  mild-slope  equations  are 


and 


^  =  0  (6.7) 

dx 


a 

dx 


(6.8) 


The  steep-slope  formulation  of  the  shallow  water  equations  are  written  in  terms  of  the 
variables  depicted  in  figure  6.33.  The  steady  one-dimensional  form  of  the  steep-slope 
equations  are 

=  0  (6,9) 

ax^ 

and 

f  ^ 

— — +  i. gd^  cos^  +  gd  sin^  =  0  (6.10) 

ax'l^  d  2  j 

where  p'  =  u'd.  For  constant  ^  the  steep-slope  formulation  is  written  in  terms  of  the 
mild-slope  formulation  variables  and  the  resulting  equation  solved  for  the  water 
surface  slope  is  given  as 
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Figure  6.33.  Schematic  of  flow  and  geometric  variables. 


dx 


Setting  C  =  0  provides  the  mild-slope  solution.  As  pointed  out  by  Berger  (1992),  the 
no  slope  case  seems  ambiguous  in  that  di^dx  =  tan  however,  di^dx  is  not  among 
the  terms  dropped  when  ^  is  assumed  to  be  small.  If  the  flow  is  supercritical,  as 
found  in  the  rectangular- to-trapezoidal  transition  test,  then  l/(cos^^  -  (Fr/cosQ^)  is  less 
than  unity  and  the  water-surface  slope  is  greater  than  the  bed  slope.  Equation  6.11 
shows  that  the  mild-slope  assumption  over-estimates  the  water-surface  slope  in  the 
case  of  supercritical  flow. 

As  the  flow  expands  and  runs  up  the  steep  side  slopes  in  the  trapezoidal 
channel,  the  mild-slope  formulation  results  in  a  water-surface  slope  that  is  too  steep 
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cos^^  - 


Fr 

cos^ 


y-J 


dzo 

dx 


(6.11) 
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which  produces  excessive  runup  on  the  side  slopes.  In  the  drawdown  case  as  the  flow 
contracts,  the  mild-slope  formulation  results  in  a  water-surface  slope  that  again  is  too 
steep  resulting  in  excessive  drawdown.  The  errors  introduced  using  the  mild-slope 
formulation  in  applications  to  supercritical  flow  are  that  on  adverse  steep  slopes  the 
flow  depth  is  too  large  and  the  velocity  too  small  and  on  favorable  steep  slopes,  the 
flow  depth  is  too  small  and  the  velocity  is  too  large.  Therefore,  the  lateral  horizontal 
distance  between  runup  and  drawdown  given  by  the  model  is  larger  than  that  observed 
in  the  flume.  Moreover,  since  the  flow  is  advected  in  the  longitudinal  direction,  this 
error  in  the  lateral  direction  (runup  and  drawdown  on  the  side  slopes)  results  in  the 
longitudinal  distance  between  wave  peaks  being  too  long.  However,  the  model  did 
capture  the  general  features  of  the  flow  and  most  importantly,  from  a  flood  control 
channel  design  standpoint,  the  maximum  wave  heights. 

Depth-averaged  longitudinal  velocity  profiles  for  three  stations  within  the  reach 
of  channel  having  sloping  sidewalls  are  provided  on  figure  6.34.  Each  plotted 
velocity,  except  those  located  near  the  waterline,  is  the  average  of  three  point 
velocities  measured  at  17,  50,  and  83  percent  of  the  flow  depth  at  a  particular  station 
and  distance  from  centerline.  At  points  near  the  waterline  a  single  measurement  of 
velocity  at  60  percent  of  the  flow  depth  was  assumed  to  represent  the  depth  average. 
Differences  between  the  model  results  and  the  flume  data  are  generally  confined  to  the 
channel  center.  The  plots  suggest  that  the  depth-averaged  velocity  of  the  flume  flow 
near  the  channel  center  is  more  uniform  across  the  channel  than  those  predicted  by  the 
model.  However,  generally  the  velocity  comparisons  are  good. 
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Figure  6.34.  Computed  and  observed  depth-averaged  velocities  for  rectangular-to- 
trapezoidal  transition  (left  and  right  directions  are  referenced  to  looking  downstream) 
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Multiple  Flow  Obstructions 

Three  flow  obstructions  composed  of  wooden  struts  were  placed  in  the 
trapezoidal  channel  having  a  base  width  of  2  ft  and  1  vertical  on  2.25  horizontal  side 
slopes  (figure  6.1),  to  simulate  the  conditions  resulting  from  supercritical  flow  through 
bridge  piers.  This  test  demonstrates  the  ease  with  which  a  finite  element  grid  can  be 
constructed  to  represent  complex  geometries  such  as  the  bridge  piers.  The  piers  were 
located  on  the  channel  base  at  the  left  and  right  toe  of  the  slopes  and  at  the  channel 
centerline.  Each  strut  was  1.0  ft  long  by  0.1  ft  wide  having  a  semicircular  nose  and 
tail  (radii  =  0.05  ft)  with  the  nose  of  each  located  at  sta  30  (figure  6.35).  A  dry  bed 
photograph  of  the  piers  installed  in  the  flume  is  provided  in  figure  6.36. 

Test  Conditions 

Tests  were  conducted  to  document  the  flow  conditions  resulting  from  a 
discharge  of  4.0  ft^/sec  with  the  tilting  flume  (figure  6.1)  set  at  a  channel  slope  of 
0.0035.  Subsequent  to  the  flume  testing,  simulations  of  these  flow  conditions  were 
conducted  to  further  evaluate  the  numerical  model’s  abilities. 

The  initial  computational  mesh  shown  on  figure  6.37a  had  1224  nodes  and 
1157  elements.  The  width  of  the  grid  was  set  based  on  the  top  width  of  flow  at 
normal  depth  in  the  absence  of  the  flow  obstructions.  Normal  depth  (0.35  ft)  was 
computed  using  a  Manning’s  roughness  coefficient  n  of  0.009  (Fr  =  1.4).  The  lateral 
mesh  resolution  at  the  upstream  end  consisted  of  2  elements  across  each  side  slope 
and  3  elements  across  the  base  of  the  trapezoidal  section.  In  the  vicinity  of  the  piers, 
the  resolution  was  increased  to  4  elements  across  each  side  slope  and  7  elements 
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Figure  6.35.  Plan  view  of  multiple  flow  obstructions  in  the  trapezoidal  channel. 

between  each  pier.  The  lower  end  of  the  model  was  represented  using  3  elements 
across  each  side  slope  and  4  elements  across  the  base.  In  a  manner  similar  to  the 
previously  reported  simulations,  the  computational  mesh  was  constructed  to  represent  a 
"headbox"  at  the  inflow  boundary  and  a  "tailbox"  at  the  outflow  boundary.  These 
nodes  and  elements  that  represented  fixed  rectangular  inflow  and  outflow  sections  and 
mild  transitions  to  and  from  trapezoidal  sections  were  added  to  the  mesh  describing 
the  trapezoidal  channel.  The  need  to  model  a  moving  boundary  representing  the 
waterline  at  the  open  boundaries  was  avoided  by  use  of  rectangular  inflow  and  outflow 
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Figure  637.  Numerical  model  computational  meshes  for  trapezoidal  channel  with 


flow  obstructions. 
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The  inflow  was  supercritical  (Froude  number  in  headbox  =  1.2)  and  therefore 
required  3  boundary  conditions.  The  grid’s  x  axis  was  aligned  with  the  channel 
centerline  so  that  the  velocity  at  the  inflow  boundary  was  in  the  x  direction.  The 
value  of  p  at  each  node  was  set  equal  to  the  unit  discharge  (total  discharge,  4.0  ft^/sec, 
divided  by  channel  width,  3.53  ft)  of  1.133  ft^/sec.  The  value  of  q  at  each  inflow 
node  was  set  to  zero.  The  inflow  depth  was  specified  as  normal  depth.  Flow  at  the 
outflow  boundary  was  supercritical  and  required  no  boundary  conditions.  The 
imposed  boundary  conditions  are  presented  in  table  6.7.  The  sidewall  boundary 
(waterline)  nodes  were  allowed  to  move  in  the  direction  of  maximum  side  slope  while 
maintaining  a  constant  depth  (fi)  of  0.01  ft  which  was  less  than  3  percent  of  normal 
depth.  The  boundaries  describing  each  pier  were  fixed,  no  flux  boundaries.  Table  6.8 
is  a  listing  of  the  model  parameters  used  in  this  simulation. 

Initial  conditions  were  specified  at  each  node  as  those  given  by  normal  depth 
flow.  The  channel  centerline  was  aligned  with  the  grid’s  x  axis  and  x-  and  y- 
directions  of  velocity  were  set  to  4.0  ft/sec  and  0.0  ft/sec,  respectively.  The  flow 
depth  at  the  channel  centerline  was  set  to  0.35  ft  such  that  the  water-surface 
longitudinal  profile  was  parallel  to  the  channel  bed  profile.  With  these  initial  and 
boundary  conditions,  simulations  were  run  using  time  steps  varying  from  0.01  sec 
(Courant  number  of  0.9)  to  0.04  sec  for  a  total  simulated  time  of  12.15  sec.  The  time 
step  was  then  increased  in  increments  to  0.2  sec  for  a  total  simulated  time  of  67.15 
sec.  At  this  time  the  solution  had  reached  steady  state. 
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Table  6.7:  Multiple  flow  obstructions,  boundary  conditions. 


Flow  Regime 

Boundary  Condition 

Supercritical 

h 

Value 


0.350  ft 


Outflow  Supercritical 


1.133  fr/sec 


0.0 


Not  Specified 


Not  Specified 


Not  Specified 


Table  6.8:  Multiple  flow  obstructions,  model  parameters. 


Condition 

Value 

a 

1.0 

p 

0.1 

n 

0.009 

0.1,  0.016  ft^/sec 

0.01  ft 
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The  grid  was  then  refined  in  the  vicinity  of  the  obstructions.  The  density  of 
the  grid  was  increased  from  immediately  upstream  of  the  piers  to  sta  36.1.  The  lateral 
resolution  on  each  of  the  side  slopes  adjacent  to  the  piers  was  increased  to  7  elements. 
The  base  of  the  trapezoidal  channel  downstream  of  the  piers  was  discretized  into  1 1 
elements.  The  resulting  computational  mesh  had  1852  nodes  and  1775  elements.  The 
initial  conditions  for  the  denser  grid  were  obtained  by  interpolation  from  the  steady 
state  solution  produced  with  the  original  grid  using  a  bilinear  interpolation  scheme. 

The  steady  state  solution  for  the  finer  mesh  was  computed  using  time  steps  varying 
from  0.02  sec  to  0.1  sec  for  an  additional  simulation  of  10.4  sec  beyond  the  steady 
state  solution  obtained  with  the  original  grid.  The  computational  mesh  associated  with 
the  steady  state  solution  is  presented  in  figure  6.37b. 

Test  Results 

Figure  6.38  shows  the  flow  conditions  produced  by  the  three  piers.  The  piers 
caused  a  choked  flow  condition  resulting  in  undulating  subcritical  flow  approaching 
the  piers.  The  flow  accelerated  between  the  obstructions  then  passed  through  critical 
depth  which  resulted  in  supercritical  flow  downstream  of  the  piers.  Figure  6.38  shows 
the  multiple  oblique  standing  waves  that  are  generated  as  the  supercritical  flow  passed 
around  the  pier  tails.  The  effects  of  these  standing  waves  on  the  shape  of  the  side 
slope  waterhnes  are  also  shown  on  the  photograph. 

Comparisons  are  made  between  the  computed  and  observed  flow  depths  by 
examination  of  the  depth  contours  presented  in  figure  6.39.  The  depths  shown  on  the 
contour  plots  are  relative  to  the  channel  centerline.  The  depth  contours  plotted  for  the 
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Figure  6.39.  Depth  contours  (feet  above  channel  bottom)  for  trapezoidal  channel  with 
flow  obstructions. 


flume  data  are  interpolations  of  a  set  of  measured  flow  depths  at  discrete  points.  The 
longitudinal  spacing  between  data  points  varied  from  as  large  as  3  ft  at  the  upstream 
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end  of  the  plotted  reach  to  0.14  ft  immediately  downstream  of  the  piers  to  0.5  ft  at  the 
downstream  end  of  the  plotted  reach.  Spacing  of  measurements  across  a  station 
ranged  from  0.5  ft  in  regions  of  smooth  water  surface  to  0.08  ft  in  the  vicinity  of  and 
just  downstream  of  the  piers.  Downstream  of  the  piers,  rapid  fluctuations  of  the  water 
surface  at  the  waterlines  presented  a  problem  in  measuring  the  time-averaged  water- 
surface  elevations  at  these  boundaries.  Although  this  error  is  not  quantified,  it  was 
greater  than  the  measurement  error  of  waterline  locations  in  the  previously  reported 
tests.  This  laboratory  error  may  explain  why  the  discrepancies  between  the  computed 
and  observed  boundaries  downstream  of  the  piers  are  larger  for  this  channel 
configuration  than  those  previously  described.  The  water-surface  profiles  provided  in 
figure  6.40  illustrate  that  the  model  and  laboratory  results  compare  quite  well  along 
the  centerline  but  that  the  waterline  in  the  flume  test  is  consistently  higher  than  the 
model  results.  This  boundary  difference  is  again  illustrated  in  the  plots  of  side 
boundaries  presented  in  figure  6.41.  Figure  6.42  is  the  computed  water-surface  mesh. 

Overall  features  of  the  flow  included  the  fact  that  although  the  channel  slope 
was  hydraulically  steep  and  would  produce  supercritical  flow  in  the  absence  of  the 
flow  obstructions,  the  obstructions  choked  the  flow.  The  specific  energy  at  normal 
flow  conditions  without  the  piers  was  less  than  that  required  to  convey  the  flow 
around  the  obstructions.  This  resulted  in  an  increase  of  specific  energy  upstream  of 
the  obstructions.  This  rise  in  depth  and  decrease  in  velocity  resulted  in  subcritical 
flow  approaching  the  obstructions.  This  phenomenon  is  commonly  referred  to  as 
Class  B  bridge  flow  (Henderson  1966). 
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STATIONS  ALONG  FLUME,  ft 


a.  Profile  along  left  wall 


b.  Profile  along  centerline 

Figure  6.40.  Water-surface  profiles  for  trapezoidal  channel  with  flow  obstructions 
(datum  is  10  ft  below  sta  0  invert  elevation). 
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Figure  6.41.  Computed  and  observed  side  boundaries  for  trapezoidal  channel  with 
flow  obstructions. 

An  undular  hydraulic  jump  formed  in  the  laboratory  flume  upstream  of  the 
obstructions.  The  flume  flow  transitioned  from  supercritical  flow  near  sta  20  to 
critical  flow  between  sta  23  and  sta  24  (critical  depth  for  this  configuration  and  flow 
rate  is  0.42  ft)  to  subcritical  flow  at  the  upstream  end  of  the  obstructions.  The  shallow 
water  equations  which  do  not  include  vertical  accelerations  are  unable  to  describe 
undular  jumps.  The  contour  plot  in  figure  6.39b  shows  that  the  numerical  model 
produced  a  hydraulic  jump  between  sta  2 1  and  sta  22.  The  contours  illustrate  a  slight 


Figure  6.42.  Computed  water-surface  mesh  for  trapezoidal  channel  with  flow 
obstructions. 


oscillation  upstream  of  the  jump.  Downstream  of  the  jump,  the  model  reproduced  the 
Sj  longitudinal  water-surface  profUe  of  the  approach  flow  as  was  observed  in  the 
flume.  An  Sj  profile  is  defined  as  decelerating  subcritical  flow  on  a  hydraulically 
steep  slope  resulting  from  downstream  control  (Henderson  1966).  Small  oscillations 
of  the  final  grid  boundary  were  noted  in  the  region  of  rapid  flow  acceleration  just 
downstream  of  the  piers  (figure  6.37b,  details  around  piers).  However,  these 
oscillations  were  confined  to  this  vicinity  and  did  not  affect  the  numerical  stability  of 
the  simulation.  The  most  important  result  from  an  engineering  standpoint  is  that  the 


153 


numerical  model  accurately  represented  the  choked  flow  condition  and  the  depth  of  the 
approaching  flow  (0.50  ft).  Accurate  predictions  of  flow  depths  are  required  for 
sidewall  height  and  minimum  bridge  soffit  elevation  determinations. 

Computed  and  observed  longitudinal  depth-averaged  velocities  are  presented  in 
figure  6.43.  Each  depth-averaged  value  plotted  for  the  flume  data  is  an  average  of 
three  velocity  measurements  taken  at  various  depths  within  the  flow  column  except  in 
the  shallow  regions  near  the  waterlines.  Velocities  were  measured  at  17,  50,  and  83 
percent  of  the  flow  depth.  Only  one  velocity  at  60  percent  of  the  depth  was  measured 
in  the  areas  of  shallow  flow. 

The  velocity  profiles  indicate  that  the  computed  flow  depth  at  sta  25.0  is  too 
large  and  the  velocities  are  too  small.  This  difference  is  also  apparent  on  the  depth 
contours  (figure  6.39).  The  model  results  indicate  that  the  flow  at  the  waterline  is  in 
the  upstream  direction.  It  was  observed  during  laboratory  testing  that  surface  currents 
near  the  banks  of  the  flume  oscillated  between  the  upstream  and  downstream 
directions,  but  the  shallow  depth  at  these  locations  prohibited  the  measurement  of 
time-averaged  velocities.  Because  the  flow  is  expanding  beginning  at  the  hydraulic 
jump,  reverse  flow  at  the  water’s  edge  of  sta  25.0  is  to  be  expected.  Near  the  end  of 
the  pier  noses  (sta  30.1)  and  downstream  of  the  piers  (sta  33.0)  the  modeled  and 
measured  velocities  are  in  close  agreement.  In  the  wake  of  the  piers  (sta  33.0)  the 
computed  longitudinal  velocities  vary  across  the  channel,  even  between  the  side  slope 
toes.  The  modeled  velocities  in  the  channel  center  vary  between  3.6  ft/sec  and  4.2 
ft/sec.  It  is  noted  that  these  velocity  variations  are  defined  by  5  nodal  values  and  are 
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LONGITUDINAL  VELOOTY,  ft/s6C 


a.  In  the  approach  flow,  sta  25.0 


b.  Near  the  pier  noses,  sta  30.1 


LONGITUDINAL  VELOCITY,  ft/sec 

c.  Downstream  of  the  piers,  sta  33.0 

Figure  6.43.  Computed  and  observed  depth-averaged  velocities  for  trapezoidal  channel 
with  flow  obstructions  (left  and  right  directions  are  referenced  to  looking  downstream). 


155 


not  node-to-node  oscillations  associated  with  numerical  instability. 

Discussion  of  Model  Results 

The  results  of  the  model  testing  demonstrate  that  the  model  can  be  used  as  a 
tool  for  the  evaluation  of  trapezoidal  high-velocity  channel  designs.  The  primary 
design  objective  of  man-made  flood  control  channels  is  the  determination  of  the  lateral 
and  longitudinal  variations  of  the  water-surface  elevation.  The  model  developed  in 
this  dissertation  study  is  sufficiently  accurate  to  serve  as  a  trapezoidal  high-velocity 
channel  design  and  evaluation  tool. 

The  test  results  also  indicate  particular  flow  features  of  trapezoidal  high- 
velocity  channels  which  the  model  cannot  reproduce.  Depth-averaged  models  (such  as 
the  one  described  in  this  study)  that  assume  uniform  vertical  distribution  of  horizontal 
velocity  are  incapable  of  reproducing  the  lateral  distribution  of  depth-averaged  velocity 
in  channel  bends  without  additional  empirical  relations  (Bernard  1993).  The  helical 
flow  in  the  real  system  results  in  the  maximum  depth-averaged  velocity  at  a  cross 
section  gradually  moving  outward  beginning  with  the  maximum  being  located  near  the 
inside  of  the  curve  at  the  curve’s  upstream  end.  Test  results  showed  that  the  depth- 
averaged  model’s  velocity  distribution  remained  constant  along  the  length  of  the  curve; 
therefore,  caution  should  be  used  in  applying  the  model  results  to  selection  of  scour 
protection  material  within  channel  bends.  Although  the  model  failed  to  accurately 
reproduce  the  velocity  distribution,  it  did  accurately  predict  the  superelevation  in  the 
water  surface  within  the  channel  bend. 

Another  significant  finding  of  the  model  testing  are  the  discrepancies  between 
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the  computed  and  observed  phases  of  the  standing  waves.  This  was  most  evident  in 
the  two  transition  cases.  In  the  trapezoidal-to-rectangular  transition  where  the  standing 
waves  were  confined  between  the  vertical  walls  of  the  rectangular  reach,  the  computed 
oblique  angle  of  the  standing  waves  was  too  large.  This  was  attributed  to  the  assumed 
hydrostatic  pressure  distribution  contained  in  the  shallow  water  equations.  However, 
downstream  of  the  rectangular-to-trapezoidal  transition  where  the  standing  waves  were 
located  within  the  trapezoidal  reach,  the  peak  values  of  the  channel  flow  width 
(distance  from  left  waterline  to  right  waterline)  trailed  those  observed  in  the  laboratory 
flume.  This  error  is  probably  due  to  a  combination  of  violations  of  the  geometrically 
mild  slope  assumption,  the  hydrostatic  pressure  distribution  assumption,  and  the 
assumption  of  no  vertical  variation  in  velocity.  The  mild-slope  formulation  of  the 
shallow  water  equations  over  predicts  runup  and  drawdown  on  the  side  slopes.  This 
overtravel  in  the  lateral  direction  also  produces  standing  wave  phase  errors  as  the  flow 
is  advected  in  the  longitudinal  direction. 

Hydrostatic  models  employing  the  mild  slope  assumption  are  not  applicable  to 
channels  having  geometrically  steep  slopes  (greater  than  about  0.1)  such  as  spillway 
chutes,  crested  drop  structures,  or  slope-faced  broad-crested  weirs.  Caution  is  also 
given  to  the  determination  of  sidewall  heights  in  or  just  downstream  of  regions 
containing  oblique  standing  waves.  The  numerical  model  results  approach  uniform 
flow  over  a  shorter  distance  than  was  observed  in  the  laboratory  flows.  Even  with  this 
short  coming,  the  model  can  be  used  to  identify  poor  designs.  The  short  transitions 
tested  in  this  study,  which  were  intentionally  designed  to  produce  large  standing 
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waves,  are  impractical  for  hydraulic  design  because  the  heights  of  the  standing  waves 
were  excessive.  The  model  could  be  used  to  evaluate  refinements  of  the  channel 
geometry  to  minimize  the  disturbance  patterns  resulting  from  the  transition. 

Positive  findings  of  the  model  testing  are  that  the  model  accurately  computed 
the  maximum  flow  depths  along  the  sidewalls  through  the  regions  of  oblique  standing 
waves  in  the  trapezoidal-to-rectangular  transition  and  the  rectangular-to-trapezoidal 
transition.  The  model  also  captured  the  magnitude  of  the  water  surface  superelevation 
in  the  horizontal  curve  of  a  trapezoidal  channel.  The  model  reproduced  the  choked 
flow  conditions  resulting  from  the  three  piers  in  the  trapezoidal  channel.  Unlike  a 
one-dimensional  analysis,  the  two-dimensional  model  does  not  require  knowledge  of 
empirical  loss  coefficients  to  describe  the  influence  of  channel  features  such  as 
transitions,  bends,  and  bridge  piers  on  the  flow  regime. 


CHAPTER  7 


SUMMARY  AND  CONCLUSIONS 
Summary 

A  two-dimensional  numerical  flow  model  for  trapezoidal  high-velocity  channels 
has  been  described.  The  model  is  designed  specifically  for  simulation  of  flow  in 
channels  having  sloping  sidewalls  in  which  the  depth  is  an  unknown  variable  in  the 
governing  equations  and  therefore,  the  plan  view  of  the  flow  domain  is  not  know  a 
priori.  Side  boundary  locations  and  initial  flow  conditions  are  assumed  and  solution  of 
the  transient  shallow  water  equations  is  advanced  in  time  until  the  steady  state  solution 
is  obtained.  As  the  computed  flow  field  evolves  to  the  steady  state,  the  model  adjusts 
the  location  of  side  boundaries  with  the  depth  solution. 

The  shallow  water  equations  modified  to  account  for  moving  boundary  effects 
are  solved  numerically  using  an  implicit  moving  finite  element  model.  The  nonlinear 
equations  are  solved  via  a  Newton-Raphson  iterative  technique  where  the  Jacobian 
consists  of  derivatives  taken  not  only  with  respect  to  the  flow  variables,  as  is  common 
in  other  flow  models,  but  also  with  respect  to  displacement  of  the  moving  boundary 
nodes.  The  test  function  used  which  is  weighted  upstream  along  characteristics  was 
chosen  because  of  its  ability  to  model  both  subcritical  and  supercritical  flow  and 
transitions  between  these  regimes. 

The  numerical  model  structure  and  solution  scheme  were  tested  against  an 
analytical  solution  of  in  viscid  shallow  water  equations  applied  to  flow  in  a  circular  V- 
shaped  channel. 

The  model  was  then  tested  by  comparing  simulation  results  with  measured 
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laboratory  data.  The  conditions  for  these  tests  were  supercritical  flows  in  a 
trapezoidal-to-rectangular  transition,  in  a  horizontal  bend  of  a  channel  having  a 
trapezoidal  cross  section,  in  a  rectangular-to-trapezoidal  transition,  and  in  a  trapezoidal 
channel  having  flow  obstructions  geometrically  similar  to  bridge  piers. 

Conclusions 

The  proposed  numerical  model  is  a  viable  means  of  simulating  flow  conditions 
in  high-velocity  channels  having  sloping  sidewalls.  This  study  produced  a  novel 
method  of  determining  the  flow  domain  and  the  flow  variables  simultaneously.  An 
alternative  approach  would  involve  updating  the  boundary  displacement  only  once  at 
the  end  of  each  time  step.  The  boundary  displacement  in  models  using  this  approach 
is  often  based  on  the  fluid  velocity  near  the  boundary  and  therefore  is  not  solved 
simultaneously  with  the  flow  depth  and  velocity.  This  method  may  be  useful  in 
subcritical  flow  where  the  nonlinear  terms  in  the  governing  equations  are  small. 
However,  in  supercritical  flow  it  is  doubtful  if  this  method  would  remain  stable  and 
convergence  of  the  solution  of  the  flow  variables  and  flow  domain  would  be  difficult 
to  achieve.  To  facilitate  convergence  and  to  avoid  instabilities,  the  boundary  location 
is  included  with  the  flow  variables  in  the  Newton-Raphson  iteration  scheme. 

Although  the  test  results  indicate  the  model’s  inability  to  reproduce  accurately 
the  phases  (wave  lengths)  of  oblique  standing  waves  in  general  and  the  lateral 
distribution  of  depth-averaged  velocity  in  a  channel  bend,  the  model  was  shown  to 
predict  accurately  the  maximum  wave  heights  occurring  in  supercritical  flow  within 
and  downstream  of  transitions,  the  superelevation  of  supercritical  flow  in  a  horizontal 
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curve  of  a  trapezoidal  channel,  and  the  choked  flow  condition  in  the  approach  to 
bridge  piers  located  in  a  trapezoidal  channel.  These  results  were  obtained  without 
empirical  coefficients,  such  as  those  used  to  account  for  geometric  features  when 
analyses  are  performed  using  one-dimensional  models.  The  model’s  ability  to 
accurately  reproduce  maximum  water  surface  elevations  is  significant  because  this  is 
the  primary  concern  to  engineers  evaluating  high-velocity  flood  control  channel 
designs. 

Although  the  two-dimensional  model  developed  and  implemented  here  is  a 
significant  advancement  in  the  numerical  modeling  of  high-velocity  channels,  it  by  no 
means  replaces  the  need  for  physical  models  of  these  hydraulic  structures.  The  model 
is  not  capable  of  reproducing  three-dimensional  flow  phenomena  nor  does  it 
incorporate  an  appropriate  turbulence  model  for  solution  of  flow  dominated  by  wall 
generated  turbulence.  Such  cases  are  found  in  channels  having  a  small  channel  width- 
to-depth  ratio  where  the  wall  boundary  friction  is  significantly  greater  than  the  bed 
friction.  Though  the  science  of  computational  fluid  dynamics  continues  rapid 
advancement,  reasons  such  as  these  lead  hydrodynamicists  and  aerodynamicists  to 
issue  the  common  cry  of  "never  abandon  the  laboratory". 

Given  these  limitations,  the  model  developed  here  can  be  most  effectively  used 
as  an  engineering  tool  to  screen  designs  prior  to  physical  laboratory  scale  modeling, 
thereby  expediting  the  design  optimization  process  and  reducing  both  the  time  and  cost 
associated  with  physical  model  testing.  Also,  the  model  can  be  beneficial  when 
applied  to  small  projects  whose  costs  do  not  justify  a  physical  model  study. 
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Recommendations  for  Additional  Research 

Development  of  schemes  for  including  the  stress  terms  in  the  shallow  water 
equations  is  an  active  current  area  of  research.  The  model  described  herein  could  be 
extended  to  include  a  better  turbulence  closure  model.  Additionally,  the  significance 
of  the  effective  stresses  arising  from  depth  averaging  of  the  equations  of  motion 
should  be  investigated,  perhaps  using  two-dimensional  momentum  correction 
coefficients.  The  bed  shear  stress  could  be  computed  using  a  Darcy-Weisbach 
formulation  for  friction  rather  than  the  Manning’s  equation  used  in  this  study.  The 
Darcy-Weisbach  friction  factor,  which  depends  on  the  Reynolds  number  and  relative 
roughness  of  the  boundary,  could  be  coded  to  account  for  variations  throughout  the 
flow  field;  however,  this  would  require  an  additional  iterative  loop  since  the  friction 
factor  and  the  flow  variables  are  interrelated.  Whether  these  model  extensions  are 
warranted  for  the  evaluation  of  flood  control  channel  designs  is  unknown. 

The  standing  wave  phase  error  could  be  reduced  using  higher  order  shallow 
water  equations.  However,  this  is  not  a  trivial  extension  and  would  necessitate  a 
reformulation  of  the  finite  element  model.  The  higher  order  partial  derivatives  found 
in  these  Boussinesq  equations  require  Hermite  polynomial  interpolation  in  which  the 
basis  function  and  its  first  derivative  are  continuous  rather  than  Lagrange  interpolation 
used  in  the  present  model  in  which  only  the  basis  function  is  continuous.  This  would 
add  to  the  computational  effort  but  computational  time  is  becoming  less  important. 

Because  the  movement  of  interior  nodes  can  be  specified  somewhat  arbitrarily, 
alternative  means  of  moving  these  nodes  could  be  addressed.  Some  sort  of  adaptive 
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regridding  scheme  could  be  incorporated,  resulting  in  highly  refined  meshes  in  areas 
where  large  gradients  of  flow  field  characteristics,  velocities  and/or  depths,  occur. 
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APPENDIX  A 


NEWTON-RAPHSON  JACOBIAN  TERMS  FOR  MOVING  FINITE  ELEMENTS 

The  Newton-Raphson  iterative  procedure  consists  of  solving  the  linear  system 
of  equations 
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for  ALj^  and  then  an  improved  estimate  for  is  obtained  from 
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Here  k  is  the  iteration  number,  j  is  the  node  location,  and 
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if  node  j  is  an  interior  node 


if  node  j  is  a  moving  boundary  node  . 
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This  procedure  is  continued  until  convergence  to  an  acceptable  residual  error  is 
obtained. 

Evaluation  of  dR^/dL^  at  a  fixed  node  consists  of  the  evaluation  of  the 
derivatives  with  respect  to  the  flow  variables  at  this  node.  This  is  commonly 
addressed  in  existing  models  and  will  not  be  discussed  here.  However,  at  a  moving 
node,  the  Jacobian  dR^/dL^^  involves  the  derivatives  dR^/ds^^,  dR^/dp^,  and 
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dR^/dq^.  An  important  aspect  of  the  model  developed  here  is  the  determination  of 
the  Jacobian  terms  for  the  finite  elements.  The  following  is  a  description  of 

these  derivatives. 

Element  variables  are  mapped  from  the  global  coordinate  system  to  local 
coordinates  (figure  4.2)  using  the  transformation 
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where  J  is  the  Jacobian  of  transformation.  Since, 
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Using  the  basis  function  (|)j,  the  derivatives  for  the  fixed  grid  case  are  interpolated  as 


where  N  is  the  number  of  nodes  on  the  element. 

The  displacement  of  node  j  is  given  as  Sj  =  (0jjj,0yj)ISjL  Then  the  global 
coordinates  are  updated  every  iteration  as 
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The  fixed  grid  derivatives  are  modified  to  account  for  nodal  displacement. 
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The  derivatives  of  the  Jacobian  of  transformation  with  respect  to  the 
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displacement  of  node  j  is 
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Gradients  of  the  basis  function,  (j),  are  given  as 
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Then  since  the  partial  derivative  of  the  element  basis  with  respect  to  global 
coordinates  is  not  a  function  of  the  displacement  of  a  global  coordinate,  i.e. 
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then  derivatives  with  respect  to  nodal  displacement  are 
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The  line  integrals  along  the  boundary  are  transformed  from  global  coordinates 
to  local  coordinates  using  the  one-dimensional  Jacobian  of  transformation,  J , 


=  i/(x2-xif  +  (yj-Yif 


(A.17) 


where  Xj  and  y^  are  the  global  coordinates  of  one  node  and  X2  and  y2  are  the  global 
coordinates  for  the  other  node.  If  node  1  needs  to  be  moved,  the  new  coordinate  is 


Xj  = 
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and  the  Jacobian  is 


176 


^  ^  (y2-yiT 

where  Ax(  =  Xj  -  Xj  -  0^^  and  AyJ  =  yj  -  yj  -  Gy^.  Then, 
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Similarly,  if  node  2  is  moved, 
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The  grid  velocity,  (Ug,  Vg),  is  approximated  as 
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where  the  superscripts  m-1,  m,  and  m+l  indicate  the  particular  time  step.  Then, 
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The  spatial  gradients  of  the  variables  are  given  as: 
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As  noted  before,  the  terms  involving  local  coordinates  are  not  a  function  of  nodal 
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displacement,  then  the  derivatives  of  the  spatial  gradients  with  respect  to  node 
movement  are 
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The  gradient  of  the  grid  velocity  is  interpolated  as 
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Then  the  derivative  with  respect  to  nodal  displacement  is 


The  length  scales  used  in  the  Petrov-Galerkin  test  function  are 
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Differentiating  with  respect  to  Sj 
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